Changes in the serum levels of IL-1 family cytokines after kidney allograft transplantation

Statistical report - data exploration and statistical analysis

Author

Cecrdlova et al.



Authors’ affiliations:

Cecrdlova E1, Krupickova L1, Fialova M1, Novotny M2, Svachova V1, Mezerova K1, Viklicky O2, Tichanek 3, Striz I1

1Department of Clinical and Transplant Immunology

2Department of Nephrology

3Department of Data Science

Institute for Clinical and Experimental Medicine, Prague, Czech Republic


The paper is under review at the Cytokine journal


1 Packages and functions

Open code
if (T) {rm(list = ls() )}
if (T) { 
  suppressWarnings(suppressMessages({
    library(tidyverse)
    library(stringr)
    library(ggpubr)
    library(emmeans)
    library(gtsummary)
    library(car)
    library(sjPlot)
    library(flextable)
    library(openxlsx)
    library(mgcv)
    library(pROC)
    library(cowplot)
    library(boot)
    library(glmnet)
    library(brms)
    library(projpred)
    library(RJDBC)
    library(janitor)
    library(arm)
    library(corrplot)
    library(lubridate)
    library(kableExtra)    
    library(ggdist)
    library(bayesplot)
    library(coxed)
    library(quantreg)
    library(rstudioapi)
    library(effects)
    library(lqmm)
    
    # Functions clashes
    select <- dplyr::select
    rename <- dplyr::rename
    mutate <- dplyr::mutate
    recode <- dplyr::recode
    summarize <- dplyr::summarize
    count <- dplyr::count
    
    # Simple math functions
    logit <-function(x){log(x/(1-x))}
    inv.logit <- function(x){exp(x)/(1+exp(x))}
  }))
}

1.1 run_model function

The function to (i) load OR (ii) run & save (if has not been done previously) results of any computation or complex table production

Open code
run <- function(expr, path, reuse = TRUE) {
  # Initialize 'fit' to NULL to ensure it's always defined
  fit <- NULL

  # Construct the file path only if 'reuse' is TRUE
  if (reuse) {
    path <- paste0(path, ".Rds")
    fit <- suppressWarnings(try(readRDS(path), silent = TRUE))

    # Check if 'fit' is an error (file not found or couldn't be read)
    if (inherits(fit, "try-error")) {
      fit <- NULL
    }
  }

  # If 'fit' is NULL (either because 'reuse' is FALSE, or the file couldn't be read), evaluate 'expr'
  if (is.null(fit)) {
    fit <- eval(substitute(expr))

    # Save 'fit' only if 'reuse' is TRUE and a valid 'path' is provided
    if (reuse && !is.null(path) && nzchar(path)) {
      saveRDS(fit, file = path)
    }
  }
  return(fit)
}

1.2 Data upload

Open code
patient_table <- read.table('data/pat_data.txt')
patient_table <- patient_table %>% select(-c(Patient, date_of_rejection,
                                             bid, Dg))
summary(patient_table)
##       id            rejection_afterTX    group           receiver_sex      
##  Length:186         Min.   :  6.00    Length:186         Length:186        
##  Class :character   1st Qu.: 17.25    Class :character   Class :character  
##  Mode  :character   Median : 69.50    Mode  :character   Mode  :character  
##                     Mean   :100.30                                         
##                     3rd Qu.:110.25                                         
##                     Max.   :505.00                                         
##                     NA's   :142                                            
##   receiver_age    donor_sex           donor_age       N_mismatch   
##  Min.   :22.00   Length:186         Min.   : 1.00   Min.   :0.000  
##  1st Qu.:45.00   Class :character   1st Qu.:44.00   1st Qu.:2.250  
##  Median :56.50   Mode  :character   Median :55.00   Median :3.000  
##  Mean   :54.67                      Mean   :52.80   Mean   :3.304  
##  3rd Qu.:65.00                      3rd Qu.:64.75   3rd Qu.:4.000  
##  Max.   :80.00                      Max.   :81.00   Max.   :6.000  
##  NA's   :48                         NA's   :48      NA's   :48     
##      dsa              anti_hla            PRA_max        PRA_actual  
##  Length:186         Length:186         Min.   : 0.00   Min.   : 0.0  
##  Class :character   Class :character   1st Qu.: 2.25   1st Qu.: 0.0  
##  Mode  :character   Mode  :character   Median :10.00   Median : 2.0  
##                                        Mean   :19.09   Mean   :11.6  
##                                        3rd Qu.:29.50   3rd Qu.:11.5  
##                                        Max.   :98.00   Max.   :92.0  
##                                        NA's   :48      NA's   :48    
##  hemodialysis_years    TX_order         ebv                cmv           
##  Min.   : 0.000     Min.   :1.000   Length:186         Length:186        
##  1st Qu.: 1.400     1st Qu.:1.000   Class :character   Class :character  
##  Median : 2.300     Median :1.000   Mode  :character   Mode  :character  
##  Mean   : 3.047     Mean   :1.109                                        
##  3rd Qu.: 4.000     3rd Qu.:1.000                                        
##  Max.   :11.900     Max.   :4.000                                        
##  NA's   :49         NA's   :48                                           
##    CIT_hours      MT_minutes    induction_therapy    tacrolimus    
##  Min.   : 0.0   Min.   : 7.00   Length:186         Min.   :0.0000  
##  1st Qu.:12.0   1st Qu.:19.00   Class :character   1st Qu.:1.0000  
##  Median :14.0   Median :24.00   Mode  :character   Median :1.0000  
##  Mean   :13.7   Mean   :25.37                      Mean   :0.9783  
##  3rd Qu.:16.0   3rd Qu.:29.75                      3rd Qu.:1.0000  
##  Max.   :24.0   Max.   :57.00                      Max.   :1.0000  
##  NA's   :48     NA's   :48                         NA's   :48      
##       MMF              mTOR           egfr_mean      creatinine_mean 
##  Min.   :0.0000   Min.   :0.00000   Min.   :0.1367   Min.   :  53.0  
##  1st Qu.:1.0000   1st Qu.:0.00000   1st Qu.:0.4558   1st Qu.: 119.0  
##  Median :1.0000   Median :0.00000   Median :0.6637   Median : 254.9  
##  Mean   :0.9275   Mean   :0.02174   Mean   :0.7720   Mean   : 277.7  
##  3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.9631   3rd Qu.: 357.2  
##  Max.   :1.0000   Max.   :1.00000   Max.   :1.8800   Max.   :1172.8  
##  NA's   :48       NA's   :48                                         
##    il1_a_mean        il1_b_mean       il1_ra_mean          il18_mean       
##  Min.   : 0.1751   Min.   : 0.3786   Min.   :    1.355   Min.   :   1.501  
##  1st Qu.: 1.4411   1st Qu.: 8.0770   1st Qu.:  667.896   1st Qu.: 278.570  
##  Median : 3.7232   Median :12.7858   Median : 1147.610   Median : 431.502  
##  Mean   : 5.0087   Mean   :12.2006   Mean   : 1624.957   Mean   : 653.351  
##  3rd Qu.: 6.5541   3rd Qu.:14.7499   3rd Qu.: 2031.337   3rd Qu.: 782.363  
##  Max.   :46.6601   Max.   :48.9390   Max.   :12548.295   Max.   :5015.730  
##                                                                            
##   il18_bp_mean   il18_free_mean     il36_b_mean     
##  Min.   : 1946   Min.   :  77.51   Min.   : 0.0258  
##  1st Qu.: 2990   1st Qu.: 206.94   1st Qu.: 2.4952  
##  Median : 3438   Median : 317.76   Median : 3.7044  
##  Mean   : 3587   Mean   : 456.12   Mean   : 5.0883  
##  3rd Qu.: 4009   3rd Qu.: 565.57   3rd Qu.: 5.7242  
##  Max.   :11804   Max.   :3279.54   Max.   :53.5479  
##  NA's   :17      NA's   :17

data_long <- read.table('data/long_data.txt')
data_long <- data_long%>% select(-Patient)
summary(data_long)
##       id             time_point           group              id_obs         
##  Length:596         Length:596         Length:596         Length:596        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##   il1_a_value       il1_b_value        il1_ra_value         il18_value      
##  Min.   : 0.1751   Min.   :  0.3786   Min.   :    1.355   Min.   :   1.501  
##  1st Qu.: 2.2933   1st Qu.: 11.1629   1st Qu.:  665.718   1st Qu.: 279.288  
##  Median : 3.7308   Median : 12.7285   Median : 1032.425   Median : 448.818  
##  Mean   : 5.5595   Mean   : 13.6417   Mean   : 1767.633   Mean   : 715.676  
##  3rd Qu.: 6.3253   3rd Qu.: 14.3496   3rd Qu.: 1745.067   3rd Qu.: 693.103  
##  Max.   :62.5197   Max.   :154.9024   Max.   :24279.375   Max.   :9421.688  
##  NA's   :96        NA's   :96         NA's   :96          NA's   :96        
##  il18_bp_value      il18_bp_ratio        il18_free        il36_b_value    
##  Min.   :   82.47   Min.   :   3.765   Min.   :  37.62   Min.   : 0.0258  
##  1st Qu.: 2976.00   1st Qu.:  31.140   1st Qu.: 193.99   1st Qu.: 2.3362  
##  Median : 3543.00   Median :  50.540   Median : 317.73   Median : 3.6714  
##  Mean   : 3708.59   Mean   :  87.717   Mean   : 486.13   Mean   : 5.5442  
##  3rd Qu.: 4168.00   3rd Qu.:  81.141   3rd Qu.: 476.33   3rd Qu.: 5.5809  
##  Max.   :30834.43   Max.   :2967.066   Max.   :6155.32   Max.   :90.0555  
##  NA's   :19         NA's   :115        NA's   :115       NA's   :96       
##    creatinine          egfr       
##  Min.   :  53.0   Min.   :0.0300  
##  1st Qu.: 114.9   1st Qu.:0.1600  
##  Median : 169.7   Median :0.5800  
##  Mean   : 315.9   Mean   :0.6238  
##  3rd Qu.: 480.1   3rd Qu.:0.9300  
##  Max.   :1778.9   Max.   :1.8800  
## 

2 Summary statistics

Summary table for all relevant variables

2.1 Healthy included

Open code
sumtab1 <- run(
 patient_table %>% select(-id, -rejection_afterTX) %>% 
  tbl_summary( by='group',
               type = list(hemodialysis_years ~ 'continuous',
                           N_mismatch ~ 'continuous')
               ) %>%
  modify_caption('Subjects characteristics according to group (healthy vs. with acute rejection vs. no acute rejection') %>% 
   add_p(),
 path = 'gitignore/run/sumtab1', reuse=TRUE)

sumtab1
Subjects characteristics according to group (healthy vs. with acute rejection vs. no acute rejection
Characteristic acute_rejection, N = 441 healthy, N = 481 no_rejection, N = 941 p-value2
receiver_sex 0.075
    female 16 (36%) 26 (54%) 33 (35%)
    male 28 (64%) 22 (46%) 61 (65%)
receiver_age 54 (44, 62) NA (NA, NA) 58 (46, 66) 0.14
    Unknown 0 48 0
donor_sex >0.9
    female 18 (41%) 0 (NA%) 37 (39%)
    male 26 (59%) 0 (NA%) 57 (61%)
    Unknown 0 48 0
donor_age 51 (42, 58) NA (NA, NA) 57 (46, 66) 0.066
    Unknown 0 48 0
N_mismatch 3 (2, 4) NA (NA, NA) 3 (3, 4) >0.9
    Unknown 0 48 0
dsa 0.002
    neg 27 (61%) 0 (NA%) 81 (86%)
    pos 17 (39%) 0 (NA%) 13 (14%)
    Unknown 0 48 0
anti_hla 0.032
    N/A 0 (0%) 0 (NA%) 1 (1.1%)
    neg 27 (61%) 0 (NA%) 75 (80%)
    pos 17 (39%) 0 (NA%) 18 (19%)
    Unknown 0 48 0
PRA_max 10 (3, 29) NA (NA, NA) 10 (2, 30) 0.9
    Unknown 0 48 0
PRA_actual 2 (0, 6) NA (NA, NA) 2 (0, 20) 0.3
    Unknown 0 48 0
hemodialysis_years 2.30 (1.40, 3.90) NA (NA, NA) 2.30 (1.30, 4.00) >0.9
    Unknown 0 48 1
TX_order <0.001
    1 32 (73%) 0 (NA%) 93 (99%)
    2 11 (25%) 0 (NA%) 1 (1.1%)
    4 1 (2.3%) 0 (NA%) 0 (0%)
    Unknown 0 48 0
ebv >0.9
    neg 0 (0%) 0 (NA%) 2 (2.8%)
    pos 35 (100%) 0 (NA%) 69 (97%)
    Unknown 9 48 23
cmv 0.2
    N/A 0 (0%) 0 (NA%) 1 (1.1%)
    neg 10 (23%) 0 (NA%) 11 (12%)
    pos 34 (77%) 0 (NA%) 82 (87%)
    Unknown 0 48 0
CIT_hours 15 (12, 18) NA (NA, NA) 14 (12, 16) 0.093
    Unknown 0 48 0
MT_minutes 22 (18, 27) NA (NA, NA) 25 (20, 31) 0.051
    Unknown 0 48 0
induction_therapy 0.004
    ATG 14 (32%) 0 (NA%) 57 (61%)
    ATG,PE,IVIG 10 (23%) 0 (NA%) 8 (8.5%)
    ATG+ Infliximab 1 (2.3%) 0 (NA%) 0 (0%)
    ATG+Infliximab 3 (6.8%) 0 (NA%) 2 (2.1%)
    Simulect 16 (36%) 0 (NA%) 27 (29%)
    Unknown 0 48 0
tacrolimus 44 (100%) 0 (NA%) 91 (97%) 0.6
    Unknown 0 48 0
MMF 39 (89%) 0 (NA%) 89 (95%) 0.3
    Unknown 0 48 0
mTOR 0 (0%) 0 (NA%) 3 (3.2%) 0.6
    Unknown 0 48 0
egfr_mean 0.49 (0.35, 0.65) 1.48 (1.15, 1.63) 0.55 (0.41, 0.72) <0.001
creatinine_mean 318 (246, 385) 85 (71, 97) 307 (235, 412) <0.001
il1_a_mean 3.6 (2.9, 4.7) 0.2 (0.2, 0.9) 5.4 (3.5, 10.2) <0.001
il1_b_mean 13.4 (12.5, 14.5) 3.8 (2.9, 4.4) 13.7 (12.2, 16.3) <0.001
il1_ra_mean 1,291 (880, 2,044) 440 (321, 577) 1,418 (1,060, 2,738) <0.001
il18_mean 626 (389, 907) 196 (146, 280) 554 (371, 1,010) <0.001
il18_bp_mean 3,389 (2,999, 3,997) 2,813 (2,465, 3,074) 3,649 (3,346, 4,194) <0.001
    Unknown 0 17 0
il18_free_mean 418 (268, 608) 145 (121, 194) 359 (236, 647) <0.001
    Unknown 0 17 0
il36_b_mean 3.5 (2.5, 4.2) 1.8 (1.3, 2.5) 5.1 (3.8, 7.2) <0.001
1 n (%); Median (IQR)
2 Pearson’s Chi-squared test; Kruskal-Wallis rank sum test; Fisher’s exact test

2.2 TX patients only

Open code
sumtab2 <- run(
 patient_table %>% select(-id, -rejection_afterTX) %>% 
   filter(group != 'healthy') %>% 
   mutate(group = factor(group)) %>% 
  tbl_summary( by='group',
               type = list(hemodialysis_years ~ 'continuous',
                           N_mismatch ~ 'continuous')
               ) %>%
  modify_caption('Patient characteristics according to developed acute rejection') %>% 
   add_p(),
 path = 'gitignore/run/sumtab2', reuse=TRUE)

sumtab2
Patient characteristics according to developed acute rejection
Characteristic acute_rejection, N = 441 no_rejection, N = 941 p-value2
receiver_sex 0.9
    female 16 (36%) 33 (35%)
    male 28 (64%) 61 (65%)
receiver_age 54 (44, 62) 58 (46, 66) 0.14
donor_sex 0.9
    female 18 (41%) 37 (39%)
    male 26 (59%) 57 (61%)
donor_age 51 (42, 58) 57 (46, 66) 0.066
N_mismatch 3 (2, 4) 3 (3, 4) >0.9
dsa <0.001
    neg 27 (61%) 81 (86%)
    pos 17 (39%) 13 (14%)
anti_hla 0.032
    N/A 0 (0%) 1 (1.1%)
    neg 27 (61%) 75 (80%)
    pos 17 (39%) 18 (19%)
PRA_max 10 (3, 29) 10 (2, 30) 0.9
PRA_actual 2 (0, 6) 2 (0, 20) 0.3
hemodialysis_years 2.30 (1.40, 3.90) 2.30 (1.30, 4.00) >0.9
    Unknown 0 1
TX_order <0.001
    1 32 (73%) 93 (99%)
    2 11 (25%) 1 (1.1%)
    4 1 (2.3%) 0 (0%)
ebv >0.9
    neg 0 (0%) 2 (2.8%)
    pos 35 (100%) 69 (97%)
    Unknown 9 23
cmv 0.2
    N/A 0 (0%) 1 (1.1%)
    neg 10 (23%) 11 (12%)
    pos 34 (77%) 82 (87%)
CIT_hours 15 (12, 18) 14 (12, 16) 0.094
MT_minutes 22 (18, 27) 25 (20, 31) 0.051
induction_therapy 0.004
    ATG 14 (32%) 57 (61%)
    ATG,PE,IVIG 10 (23%) 8 (8.5%)
    ATG+ Infliximab 1 (2.3%) 0 (0%)
    ATG+Infliximab 3 (6.8%) 2 (2.1%)
    Simulect 16 (36%) 27 (29%)
tacrolimus 44 (100%) 91 (97%) 0.6
MMF 39 (89%) 89 (95%) 0.3
mTOR 0 (0%) 3 (3.2%) 0.6
egfr_mean 0.49 (0.35, 0.65) 0.55 (0.41, 0.72) 0.2
creatinine_mean 318 (246, 385) 307 (235, 412) 0.9
il1_a_mean 3.6 (2.9, 4.7) 5.4 (3.5, 10.2) 0.001
il1_b_mean 13.4 (12.5, 14.5) 13.7 (12.2, 16.3) 0.13
il1_ra_mean 1,291 (880, 2,044) 1,418 (1,060, 2,738) 0.10
il18_mean 626 (389, 907) 554 (371, 1,010) >0.9
il18_bp_mean 3,389 (2,999, 3,997) 3,649 (3,346, 4,194) 0.2
il18_free_mean 418 (268, 608) 359 (236, 647) 0.7
il36_b_mean 3.5 (2.5, 4.2) 5.1 (3.8, 7.2) <0.001
1 n (%); Median (IQR)
2 Pearson’s Chi-squared test; Wilcoxon rank sum test; Fisher’s exact test

2.3 Summary table for continuous variables

Open code
sumtab <- patient_table %>%
  select(
    receiver_age, donor_age,
    N_mismatch, PRA_max, PRA_actual, CIT_hours, MT_minutes,
    hemodialysis_years, group, il1_a_mean, il1_b_mean, il1_ra_mean,
    il18_mean, il18_bp_mean, il18_free_mean, il36_b_mean) %>% 
  group_by(group) %>%
  summarize(across(where(is.numeric), 
                   list(median = ~median(., na.rm = TRUE), 
                        mean = ~mean(., na.rm = TRUE), 
                        IQR = ~IQR(., na.rm = TRUE), 
                        Q25 = ~quantile(., 0.25, na.rm = TRUE), 
                        Q75 = ~quantile(., 0.75, na.rm = TRUE),
                        sd = ~sd(., na.rm = TRUE)), 
                   .names = "{.col}_{.fn}")) 
kableExtra::kable(t(sumtab)) 
group acute_rejection healthy no_rejection
receiver_age_median 54 NA 58
receiver_age_mean 52.15909 NA 55.85106
receiver_age_IQR 18.25 NA 19.75
receiver_age_Q25 44.00 NA 46.25
receiver_age_Q75 62.25 NA 66.00
receiver_age_sd 13.95498 NA 13.41076
donor_age_median 51.0 NA 56.5
donor_age_mean 49.90909 NA 54.15957
donor_age_IQR 16.5 NA 20.0
donor_age_Q25 41.75 NA 46.00
donor_age_Q75 58.25 NA 66.00
donor_age_sd 14.00468 NA 15.41630
N_mismatch_median 3 NA 3
N_mismatch_mean 3.250000 NA 3.329787
N_mismatch_IQR 2 NA 1
N_mismatch_Q25 2 NA 3
N_mismatch_Q75 4 NA 4
N_mismatch_sd 1.331637 NA 1.212731
PRA_max_median 10 NA 10
PRA_max_mean 19.40909 NA 18.93617
PRA_max_IQR 26.25 NA 27.25
PRA_max_Q25 2.75 NA 2.25
PRA_max_Q75 29.0 NA 29.5
PRA_max_sd 23.26322 NA 19.97407
PRA_actual_median 2 NA 2
PRA_actual_mean 10.00000 NA 12.35106
PRA_actual_IQR 6.0 NA 19.5
PRA_actual_Q25 0 NA 0
PRA_actual_Q75 6.0 NA 19.5
PRA_actual_sd 21.58703 NA 18.72191
CIT_hours_median 15 NA 14
CIT_hours_mean 14.31818 NA 13.40426
CIT_hours_IQR 6 NA 4
CIT_hours_Q25 12 NA 12
CIT_hours_Q75 18 NA 16
CIT_hours_sd 5.065950 NA 4.489699
MT_minutes_median 21.5 NA 25.0
MT_minutes_mean 23.90909 NA 26.05319
MT_minutes_IQR 8.50 NA 10.75
MT_minutes_Q25 18.00 NA 20.25
MT_minutes_Q75 26.5 NA 31.0
MT_minutes_sd 8.454580 NA 8.773576
hemodialysis_years_median 2.3 NA 2.3
hemodialysis_years_mean 2.972727 NA 3.082796
hemodialysis_years_IQR 2.5 NA 2.7
hemodialysis_years_Q25 1.4 NA 1.3
hemodialysis_years_Q75 3.9 NA 4.0
hemodialysis_years_sd 2.271396 NA 2.491273
il1_a_mean_median 3.6138431 0.1751254 5.4099081
il1_a_mean_mean 4.519618 1.790396 6.880929
il1_a_mean_IQR 1.8399964 0.7684427 6.6405988
il1_a_mean_Q25 2.8678507 0.1751254 3.5186279
il1_a_mean_Q75 4.7078471 0.9435681 10.1592267
il1_a_mean_sd 3.702535 6.758296 5.079430
il1_b_mean_median 13.379882 3.800578 13.726654
il1_b_mean_mean 13.614621 4.568534 15.435899
il1_b_mean_IQR 2.016068 1.490545 4.087478
il1_b_mean_Q25 12.456512 2.865737 12.211393
il1_b_mean_Q75 14.472580 4.356282 16.298871
il1_b_mean_sd 1.876318 6.136904 5.552799
il1_ra_mean_median 1291.0082 439.8988 1417.6213
il1_ra_mean_mean 1607.137 531.280 2191.771
il1_ra_mean_IQR 1163.9840 255.8475 1677.9097
il1_ra_mean_Q25 879.6144 320.7617 1059.8544
il1_ra_mean_Q75 2043.5983 576.6092 2737.7641
il1_ra_mean_sd 1054.4684 464.7765 1883.0248
il18_mean_median 626.2038 196.2464 553.9956
il18_mean_mean 706.7674 220.7274 849.2618
il18_mean_IQR 518.0835 133.1681 639.2590
il18_mean_Q25 388.6040 146.4606 370.9122
il18_mean_Q75 906.6875 279.6288 1010.1712
il18_mean_sd 413.15414 99.47696 778.91956
il18_bp_mean_median 3389.167 2813.000 3649.250
il18_bp_mean_mean 3666.610 2854.484 3791.695
il18_bp_mean_IQR 998.5208 609.0000 847.5291
il18_bp_mean_Q25 2998.875 2465.000 3346.188
il18_bp_mean_Q75 3997.396 3074.000 4193.717
il18_bp_mean_sd 889.4910 514.8316 1053.8108
il18_free_mean_median 417.7738 144.8761 358.7065
il18_free_mean_mean 478.2362 164.6453 541.8925
il18_free_mean_IQR 339.20872 72.47453 411.11565
il18_free_mean_Q25 268.4762 121.4914 235.7412
il18_free_mean_Q75 607.6849 193.9659 646.8569
il18_free_mean_sd 286.58937 69.54324 499.73260
il36_b_mean_median 3.516403 1.823203 5.079123
il36_b_mean_mean 5.616041 2.465567 6.180597
il36_b_mean_IQR 1.639928 1.250801 3.346720
il36_b_mean_Q25 2.543047 1.279165 3.806194
il36_b_mean_Q75 4.182975 2.529966 7.152914
il36_b_mean_sd 8.701474 3.398531 4.137029

3 Cytokine levels changes over time and across groups

Creating a new variable, group_sep, that combines group and time_point into a single variable. Re-defining group and time_point to be factors

Open code
data_long <- data_long %>% 
  mutate(group_sep = factor(interaction(group, time_point)),
         group = factor(group, levels = 
                          c("healthy",
                            "no_rejection", 
                            "acute_rejection")),
         time_point = factor(time_point))

3.1 Il1_a

3.1.1 Mixed model

Open code

## filtering data
model_data <- data_long %>% 
  filter(!is.na(il1_a_value))

## model fit
model_il1_a <-  lmer(log10(il1_a_value) ~ 
                       group_sep +
                       (1|id), data = model_data)

## diagnostic plot of model
plot(model_il1_a)

Open code

## anova result od model
Anova(model_il1_a, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_a_value)
##                F Df Df.res    Pr(>F)    
## group_sep 15.739  9 417.45 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## model summary
summary(model_il1_a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_a_value) ~ group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 807.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5688 -0.3288  0.0750  0.5766  3.4928 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.05345  0.2312  
##  Residual             0.23773  0.4876  
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      0.54459    0.08602   6.331
## group_sephealthy.000            -0.96152    0.11604  -8.286
## group_sepno_rejection.000        0.01385    0.10245   0.135
## group_sepacute_rejection.007    -0.08337    0.10815  -0.771
## group_sepno_rejection.007       -0.09297    0.10245  -0.907
## group_sepacute_rejection.090    -0.02344    0.10975  -0.214
## group_sepno_rejection.090        0.19470    0.12090   1.610
## group_sepacute_rejection.365    -0.11856    0.11488  -1.032
## group_sepno_rejection.365        0.13162    0.12090   1.089
## group_sepacute_rejection.biopsy -0.06165    0.12267  -0.503
## 
## Correlation of Fixed Effects:
##              (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000  -0.741                                                      
## grp_sp_.000  -0.840  0.622                                               
## grp_spc_.007 -0.663  0.491  0.557                                        
## grp_spn_.007 -0.840  0.622  0.759   0.557                                
## grp_spc_.090 -0.656  0.486  0.551   0.520        0.551                   
## grp_spn_.090 -0.711  0.527  0.643   0.472        0.643        0.467      
## grp_spc_.365 -0.628  0.466  0.528   0.498        0.528        0.492      
## grp_spn_.365 -0.711  0.527  0.643   0.472        0.643        0.467      
## grp_spct_r.  -0.586  0.434  0.492   0.464        0.492        0.459      
##              grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000                                        
## grp_sp_.000                                        
## grp_spc_.007                                       
## grp_spn_.007                                       
## grp_spc_.090                                       
## grp_spn_.090                                       
## grp_spc_.365  0.447                                
## grp_spn_.365  0.583        0.447                   
## grp_spct_r.   0.417        0.435        0.417

3.1.2 Pairwise comparison

3.1.2.1 List of contrasts

List of contrast names and left and right coordinates on the X axis

Open code
contrast <- c(

  ## healthy vs. others
  "healthy vs. no rejection D_0",
  "healthy vs. no rejection D_7",
  "healthy vs. no rejection D_90",
  "healthy vs. no rejection D_365",
  "healthy vs. acute rejection D_0",
  "healthy vs. acute rejection D_7",
  "healthy vs. acute rejection D_90",
  "healthy vs. acute rejection D_365",
  "healthy vs. acute rejection D_biopsy",


  ## between-subjects comparisons
  "no rejection D_0 vs. acute rejection D_0",
  "no rejection D_7 vs. acute rejection D_7",
  "no rejection D_90 vs. acute rejection D_90",
  "no rejection D_365 vs. acute rejection D_365",
  "no rejection D_0 vs. acute rejection D_biopsy",
  "no rejection D_7 vs. acute rejection D_biopsy",
  "no rejection D_90 vs. acute rejection D_biopsy",
  "no rejection D_365 vs. acute rejection D_biopsy",

  ## within-subjects comparisons
  ### no rejection
  "no rejection D_0 vs. no rejection D_7",
  "no rejection D_0 vs. no rejection D_90",
  "no rejection D_0 vs. no rejection D_365",
  "no rejection D_7 vs. no rejection D_90",
  "no rejection D_7 vs. no rejection D_365",
  "no rejection D_90 vs. no rejection D_365",

    ### acute rejection
  "acute rejection D_0 vs. acute rejection D_7",
  "acute rejection D_0 vs. acute rejection D_90",
  "acute rejection D_0 vs. acute rejection D_365",
  "acute rejection D_7 vs. acute rejection D_90",
  "acute rejection D_7 vs. acute rejection D_365",
  "acute rejection D_90 vs. acute rejection D_365",
  "acute rejection D_0 vs. acute rejection D_biopsy",
  "acute rejection D_7 vs. acute rejection D_biopsy",
  "acute rejection D_90 vs. acute rejection D_biopsy",
  "acute rejection D_365 vs. acute rejection D_biopsy"
  )
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
              rep(0.8, 3), rep(1.8, 2), 2.8,
              rep(1.2, 3), rep(2.2, 2), 3.2,
              seq(1.2, 4.2, 1))

xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
              5, seq(1.2, 4.2, 1), rep(5, 4),
              seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
              2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
              rep(5, 4))

anot_coord <- (xcoord_1+xcoord_2)/2

pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
                         anot_coord)

pair_table <- pair_table %>% 
  mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>% 
  mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
  mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>% 
  mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))

  
  pair_table$p <- NA

3.1.2.2 Between-subjects comparisons

Open code
set.seed(16)
for (i in 1:17){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i])|
           group_sep == paste(pair_table$term2[i]))
  
  g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
  g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
  
  pair_table$p[i] <- wilcox.test(g1$il1_a_value, g2$il1_a_value, paired=FALSE)$p.value
}

3.1.2.3 Within-subject comparisons

Open code
set.seed(16)
for (i in 18:33){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i]) |
           group_sep == paste(pair_table$term2[i])) %>% 
    select(id, group_sep, il1_a_value) %>% 
    pivot_wider(names_from = group_sep, values_from = il1_a_value) %>% 
    filter(complete.cases(.)) %>% data.frame()
  
  g1 <- dat_pair[,2]
  g2 <- dat_pair[,3]
  
  pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}

pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),] %>% data.frame()
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')

3.1.2.4 Resulting table

Open code
## the whole `pair table`
pair_table
##                                           contrast xcoord_1 xcoord_2 anot_coord
## 1                     healthy vs. no rejection D_0      6.0      0.8        3.4
## 2                     healthy vs. no rejection D_7      6.0      1.8        3.9
## 3                    healthy vs. no rejection D_90      6.0      2.8        4.4
## 4                   healthy vs. no rejection D_365      6.0      3.8        4.9
## 5                  healthy vs. acute rejection D_0      6.0      1.2        3.6
## 6                  healthy vs. acute rejection D_7      6.0      2.2        4.1
## 7                 healthy vs. acute rejection D_90      6.0      3.2        4.6
## 8                healthy vs. acute rejection D_365      6.0      4.2        5.1
## 9             healthy vs. acute rejection D_biopsy      6.0      5.0        5.5
## 10        no rejection D_0 vs. acute rejection D_0      0.8      1.2        1.0
## 11        no rejection D_7 vs. acute rejection D_7      1.8      2.2        2.0
## 12      no rejection D_90 vs. acute rejection D_90      2.8      3.2        3.0
## 13    no rejection D_365 vs. acute rejection D_365      3.8      4.2        4.0
## 16  no rejection D_90 vs. acute rejection D_biopsy      2.8      5.0        3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy      3.8      5.0        4.4
## 18           no rejection D_0 vs. no rejection D_7      0.8      1.8        1.3
## 21          no rejection D_7 vs. no rejection D_90      1.8      2.8        2.3
## 23        no rejection D_90 vs. no rejection D_365      2.8      3.8        3.3
## 24     acute rejection D_0 vs. acute rejection D_7      1.2      2.2        1.7
## 27    acute rejection D_7 vs. acute rejection D_90      2.2      3.2        2.7
## 29  acute rejection D_90 vs. acute rejection D_365      3.2      4.2        3.7
##                  term1                  term2            p  cole          fdr
## 1          healthy.000       no_rejection.000 6.973501e-11 black 2.928870e-10
## 2          healthy.000       no_rejection.007 5.227359e-11 black 2.744363e-10
## 3          healthy.000       no_rejection.090 1.475037e-13 black 3.097579e-12
## 4          healthy.000       no_rejection.365 7.054097e-13 black 7.406802e-12
## 5          healthy.000    acute_rejection.000 5.467856e-10 black 1.913750e-09
## 6          healthy.000    acute_rejection.007 2.952645e-09 black 7.750693e-09
## 7          healthy.000    acute_rejection.090 3.797214e-11 black 2.658050e-10
## 8          healthy.000    acute_rejection.365 2.847456e-09 black 7.750693e-09
## 9          healthy.000 acute_rejection.biopsy 4.447710e-09 black 1.037799e-08
## 10    no_rejection.000    acute_rejection.000 3.608654e-01 black 4.300071e-01
## 11    no_rejection.007    acute_rejection.007 3.455284e-01 black 4.300071e-01
## 12    no_rejection.090    acute_rejection.090 6.924033e-06 black 1.321861e-05
## 13    no_rejection.365    acute_rejection.365 2.968675e-06 black 6.234218e-06
## 16    no_rejection.090 acute_rejection.biopsy 1.026038e-04 black 1.795567e-04
## 17    no_rejection.365 acute_rejection.biopsy 7.185841e-04 black 1.160790e-03
## 18    no_rejection.000       no_rejection.007 7.523005e-03 blue3 1.128451e-02
## 21    no_rejection.007       no_rejection.090 6.940287e-01 blue3 7.670843e-01
## 23    no_rejection.090       no_rejection.365 3.685775e-01 blue3 4.300071e-01
## 24 acute_rejection.000    acute_rejection.007 8.076073e-01  red4 8.076073e-01
## 27 acute_rejection.007    acute_rejection.090 7.446044e-01  red4 7.818346e-01
## 29 acute_rejection.090    acute_rejection.365 1.794364e-01  red4 2.512110e-01

## only significant differences
pair_table %>% 
  filter(p<0.05) %>% 
  select(contrast, p, fdr) %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()
##                                           contrast            p          fdr
## 1                     healthy vs. no rejection D_0 6.973501e-11 2.928870e-10
## 2                     healthy vs. no rejection D_7 5.227359e-11 2.744363e-10
## 3                    healthy vs. no rejection D_90 1.475037e-13 3.097579e-12
## 4                   healthy vs. no rejection D_365 7.054097e-13 7.406802e-12
## 5                  healthy vs. acute rejection D_0 5.467856e-10 1.913750e-09
## 6                  healthy vs. acute rejection D_7 2.952645e-09 7.750693e-09
## 7                 healthy vs. acute rejection D_90 3.797214e-11 2.658050e-10
## 8                healthy vs. acute rejection D_365 2.847456e-09 7.750693e-09
## 9             healthy vs. acute rejection D_biopsy 4.447710e-09 1.037799e-08
## 10      no rejection D_90 vs. acute rejection D_90 6.924033e-06 1.321861e-05
## 11    no rejection D_365 vs. acute rejection D_365 2.968675e-06 6.234218e-06
## 12  no rejection D_90 vs. acute rejection D_biopsy 1.026038e-04 1.795567e-04
## 13 no rejection D_365 vs. acute rejection D_biopsy 7.185841e-04 1.160790e-03
## 14           no rejection D_0 vs. no rejection D_7 7.523005e-03 1.128451e-02
##    star
## 1   ***
## 2   ***
## 3   ***
## 4   ***
## 5   ***
## 6   ***
## 7   ***
## 8   ***
## 9   ***
## 10  ***
## 11  ***
## 12  ***
## 13  ***
## 14   **


pair_rest <- pair_table %>% 
  filter(p<0.05,
         term1 != 'healthy.000') %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()

3.1.3 Visualization

Open code

## plot
plotos <- data_long %>% 
  filter(!is.na(il1_a_value)) %>% 
  mutate(group = factor(group, levels = 
                          c("no_rejection", 
                            "acute_rejection",
                            "healthy")),
         time_point = factor(time_point)) %>% 
  mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>% 
  mutate(group_sep = interaction(group, time_point)) %>% 
  
## plotting
  ggplot(aes(x=time_point, y=il1_a_value)) +
  geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA, 
               size=0.5,position = position_dodge(width = 0.84), color='black') +
  labs(x = "Time post TX", y = expression("IL-1" * alpha * " (pg/ml)")) +
  scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
                              "Biopsy", "Healthy")) +
  theme(axis.text = element_text(size = 11),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 12)) +
  scale_fill_manual(values = c("no_rejection" = "skyblue2",
                               "acute_rejection" = "coral",
                               "healthy" = "darkolivegreen3"), 
                    labels = c("no rejection", "acute rejection", "healthy")) +
    coord_cartesian(ylim = c(0, 30)) +
  annotate("text", x = 6, y = 5, label = "#", size = 7) 

pair_rest$y <- c(12, 11, 16, 13,30)

plotos <- plotos + geom_segment(data = pair_rest, 
                                aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
                                color = pair_rest$cole,
                                linetype = "solid", size = 0.4) 
plot_il1_a <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+0.5, 
                            label = pair_rest$star, size = 7, color=pair_rest$cole) 
plot_il1_a

3.2 Il1_b

3.2.1 Mixed model

Open code

## filtering data
model_data <- data_long %>% 
  filter(!is.na(il1_b_value))

## model fit
model_il1_b <-  lmer(log10(il1_b_value) ~ 
                       group_sep +
                       (1|id), data = model_data)

## diagnostic plot of model
plot(model_il1_b)

Open code

## anova result od model
Anova(model_il1_b, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_b_value)
##               F Df Df.res    Pr(>F)    
## group_sep 83.17  9 417.24 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## model summary
summary(model_il1_b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_b_value) ~ group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: -509.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.2908 -0.5140 -0.0895  0.2629  7.1938 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.003725 0.06103 
##  Residual             0.016117 0.12695 
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      1.16016    0.02245  51.675
## group_sephealthy.000            -0.59874    0.03029 -19.768
## group_sepno_rejection.000        0.04741    0.02674   1.773
## group_sepacute_rejection.007    -0.05044    0.02816  -1.791
## group_sepno_rejection.007       -0.01361    0.02674  -0.509
## group_sepacute_rejection.090    -0.04086    0.02858  -1.430
## group_sepno_rejection.090       -0.07554    0.03155  -2.395
## group_sepacute_rejection.365    -0.05930    0.02992  -1.982
## group_sepno_rejection.365       -0.09149    0.03155  -2.900
## group_sepacute_rejection.biopsy -0.04640    0.03195  -1.453
## 
## Correlation of Fixed Effects:
##              (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000  -0.741                                                      
## grp_sp_.000  -0.840  0.622                                               
## grp_spc_.007 -0.661  0.490  0.555                                        
## grp_spn_.007 -0.840  0.622  0.760   0.555                                
## grp_spc_.090 -0.654  0.485  0.549   0.520        0.549                   
## grp_spn_.090 -0.712  0.528  0.644   0.471        0.644        0.466      
## grp_spc_.365 -0.627  0.465  0.526   0.498        0.526        0.492      
## grp_spn_.365 -0.712  0.528  0.644   0.471        0.644        0.466      
## grp_spct_r.  -0.584  0.433  0.490   0.464        0.490        0.459      
##              grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000                                        
## grp_sp_.000                                        
## grp_spc_.007                                       
## grp_spn_.007                                       
## grp_spc_.090                                       
## grp_spn_.090                                       
## grp_spc_.365  0.446                                
## grp_spn_.365  0.585        0.446                   
## grp_spct_r.   0.416        0.435        0.416

3.2.2 Pairwise comparison

3.2.2.1 List of contrasts

List of contrast names and left and right coordinates on the X axis

Open code
contrast <- c(

  ## healthy vs. others
  "healthy vs. no rejection D_0",
  "healthy vs. no rejection D_7",
  "healthy vs. no rejection D_90",
  "healthy vs. no rejection D_365",
  "healthy vs. acute rejection D_0",
  "healthy vs. acute rejection D_7",
  "healthy vs. acute rejection D_90",
  "healthy vs. acute rejection D_365",
  "healthy vs. acute rejection D_biopsy",


  ## between-subjects comparisons
  "no rejection D_0 vs. acute rejection D_0",
  "no rejection D_7 vs. acute rejection D_7",
  "no rejection D_90 vs. acute rejection D_90",
  "no rejection D_365 vs. acute rejection D_365",
  "no rejection D_0 vs. acute rejection D_biopsy",
  "no rejection D_7 vs. acute rejection D_biopsy",
  "no rejection D_90 vs. acute rejection D_biopsy",
  "no rejection D_365 vs. acute rejection D_biopsy",

  ## within-subjects comparisons
  ### no rejection
  "no rejection D_0 vs. no rejection D_7",
  "no rejection D_0 vs. no rejection D_90",
  "no rejection D_0 vs. no rejection D_365",
  "no rejection D_7 vs. no rejection D_90",
  "no rejection D_7 vs. no rejection D_365",
  "no rejection D_90 vs. no rejection D_365",

    ### acute rejection
  "acute rejection D_0 vs. acute rejection D_7",
  "acute rejection D_0 vs. acute rejection D_90",
  "acute rejection D_0 vs. acute rejection D_365",
  "acute rejection D_7 vs. acute rejection D_90",
  "acute rejection D_7 vs. acute rejection D_365",
  "acute rejection D_90 vs. acute rejection D_365",
  "acute rejection D_0 vs. acute rejection D_biopsy",
  "acute rejection D_7 vs. acute rejection D_biopsy",
  "acute rejection D_90 vs. acute rejection D_biopsy",
  "acute rejection D_365 vs. acute rejection D_biopsy"
  )
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
              rep(0.8, 3), rep(1.8, 2), 2.8,
              rep(1.2, 3), rep(2.2, 2), 3.2,
              seq(1.2, 4.2, 1))

xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
              5, seq(1.2, 4.2, 1), rep(5, 4),
              seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
              2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
              rep(5, 4))

anot_coord <- (xcoord_1+xcoord_2)/2

pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
                         anot_coord)

pair_table <- pair_table %>% 
  mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>% 
  mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
  mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>% 
  mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
  
  pair_table$p <- NA

3.2.2.2 Between-subjects comparisons

Open code
set.seed(16)
for (i in 1:17){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i])|
           group_sep == paste(pair_table$term2[i]))
  
  g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
  g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
  
  pair_table$p[i] <- wilcox.test(g1$il1_b_value, g2$il1_b_value, paired=FALSE)$p.value
}

3.2.2.3 Within-subject comparisons

Open code
set.seed(16)
for (i in 18:33){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i]) |
           group_sep == paste(pair_table$term2[i])) %>% 
    select(id, group_sep, il1_b_value) %>% 
    pivot_wider(names_from = group_sep, values_from = il1_b_value) %>% 
    filter(complete.cases(.)) %>% data.frame()
  
  g1 <- dat_pair[,2]
  g2 <- dat_pair[,3]
  
  pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}


pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')

3.2.2.4 Resulting table

Open code
## the whole `pair table`
pair_table
##                                           contrast xcoord_1 xcoord_2 anot_coord
## 1                     healthy vs. no rejection D_0      6.0      0.8        3.4
## 2                     healthy vs. no rejection D_7      6.0      1.8        3.9
## 3                    healthy vs. no rejection D_90      6.0      2.8        4.4
## 4                   healthy vs. no rejection D_365      6.0      3.8        4.9
## 5                  healthy vs. acute rejection D_0      6.0      1.2        3.6
## 6                  healthy vs. acute rejection D_7      6.0      2.2        4.1
## 7                 healthy vs. acute rejection D_90      6.0      3.2        4.6
## 8                healthy vs. acute rejection D_365      6.0      4.2        5.1
## 9             healthy vs. acute rejection D_biopsy      6.0      5.0        5.5
## 10        no rejection D_0 vs. acute rejection D_0      0.8      1.2        1.0
## 11        no rejection D_7 vs. acute rejection D_7      1.8      2.2        2.0
## 12      no rejection D_90 vs. acute rejection D_90      2.8      3.2        3.0
## 13    no rejection D_365 vs. acute rejection D_365      3.8      4.2        4.0
## 16  no rejection D_90 vs. acute rejection D_biopsy      2.8      5.0        3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy      3.8      5.0        4.4
## 18           no rejection D_0 vs. no rejection D_7      0.8      1.8        1.3
## 21          no rejection D_7 vs. no rejection D_90      1.8      2.8        2.3
## 23        no rejection D_90 vs. no rejection D_365      2.8      3.8        3.3
## 24     acute rejection D_0 vs. acute rejection D_7      1.2      2.2        1.7
## 27    acute rejection D_7 vs. acute rejection D_90      2.2      3.2        2.7
## 29  acute rejection D_90 vs. acute rejection D_365      3.2      4.2        3.7
##                  term1                  term2            p  cole          fdr
## 1          healthy.000       no_rejection.000 1.007494e-20 black 1.100496e-19
## 2          healthy.000       no_rejection.007 1.048092e-20 black 1.100496e-19
## 3          healthy.000       no_rejection.090 1.941165e-14 black 6.224985e-14
## 4          healthy.000       no_rejection.365 2.074995e-14 black 6.224985e-14
## 5          healthy.000    acute_rejection.000 1.954387e-14 black 6.224985e-14
## 6          healthy.000    acute_rejection.007 3.839545e-15 black 2.687682e-14
## 7          healthy.000    acute_rejection.090 8.457201e-15 black 4.440030e-14
## 8          healthy.000    acute_rejection.365 1.162541e-13 black 3.051669e-13
## 9          healthy.000 acute_rejection.biopsy 4.140081e-12 black 9.660189e-12
## 10    no_rejection.000    acute_rejection.000 1.656273e-01 black 2.045985e-01
## 11    no_rejection.007    acute_rejection.007 3.029547e-02 black 4.241366e-02
## 12    no_rejection.090    acute_rejection.090 1.269805e-01 black 1.666619e-01
## 13    no_rejection.365    acute_rejection.365 1.059940e-02 black 1.712211e-02
## 16    no_rejection.090 acute_rejection.biopsy 2.524270e-01 black 2.944981e-01
## 17    no_rejection.365 acute_rejection.biopsy 2.029089e-02 black 3.043634e-02
## 18    no_rejection.000       no_rejection.007 8.577426e-05 blue3 1.801260e-04
## 21    no_rejection.007       no_rejection.090 1.573119e-03 blue3 3.003227e-03
## 23    no_rejection.090       no_rejection.365 5.391972e-01 blue3 5.391972e-01
## 24 acute_rejection.000    acute_rejection.007 5.187312e-03  red4 9.077796e-03
## 27 acute_rejection.007    acute_rejection.090 5.013747e-01  red4 5.264435e-01
## 29 acute_rejection.090    acute_rejection.365 4.890153e-01  red4 5.264435e-01

## only significant differences
pair_table %>% filter(p<0.05) %>% data.frame()
##                                           contrast xcoord_1 xcoord_2 anot_coord
## 1                     healthy vs. no rejection D_0      6.0      0.8        3.4
## 2                     healthy vs. no rejection D_7      6.0      1.8        3.9
## 3                    healthy vs. no rejection D_90      6.0      2.8        4.4
## 4                   healthy vs. no rejection D_365      6.0      3.8        4.9
## 5                  healthy vs. acute rejection D_0      6.0      1.2        3.6
## 6                  healthy vs. acute rejection D_7      6.0      2.2        4.1
## 7                 healthy vs. acute rejection D_90      6.0      3.2        4.6
## 8                healthy vs. acute rejection D_365      6.0      4.2        5.1
## 9             healthy vs. acute rejection D_biopsy      6.0      5.0        5.5
## 10        no rejection D_7 vs. acute rejection D_7      1.8      2.2        2.0
## 11    no rejection D_365 vs. acute rejection D_365      3.8      4.2        4.0
## 12 no rejection D_365 vs. acute rejection D_biopsy      3.8      5.0        4.4
## 13           no rejection D_0 vs. no rejection D_7      0.8      1.8        1.3
## 14          no rejection D_7 vs. no rejection D_90      1.8      2.8        2.3
## 15     acute rejection D_0 vs. acute rejection D_7      1.2      2.2        1.7
##                  term1                  term2            p  cole          fdr
## 1          healthy.000       no_rejection.000 1.007494e-20 black 1.100496e-19
## 2          healthy.000       no_rejection.007 1.048092e-20 black 1.100496e-19
## 3          healthy.000       no_rejection.090 1.941165e-14 black 6.224985e-14
## 4          healthy.000       no_rejection.365 2.074995e-14 black 6.224985e-14
## 5          healthy.000    acute_rejection.000 1.954387e-14 black 6.224985e-14
## 6          healthy.000    acute_rejection.007 3.839545e-15 black 2.687682e-14
## 7          healthy.000    acute_rejection.090 8.457201e-15 black 4.440030e-14
## 8          healthy.000    acute_rejection.365 1.162541e-13 black 3.051669e-13
## 9          healthy.000 acute_rejection.biopsy 4.140081e-12 black 9.660189e-12
## 10    no_rejection.007    acute_rejection.007 3.029547e-02 black 4.241366e-02
## 11    no_rejection.365    acute_rejection.365 1.059940e-02 black 1.712211e-02
## 12    no_rejection.365 acute_rejection.biopsy 2.029089e-02 black 3.043634e-02
## 13    no_rejection.000       no_rejection.007 8.577426e-05 blue3 1.801260e-04
## 14    no_rejection.007       no_rejection.090 1.573119e-03 blue3 3.003227e-03
## 15 acute_rejection.000    acute_rejection.007 5.187312e-03  red4 9.077796e-03

pair_rest <- pair_table %>% 
  filter(p<0.05,
         term1 != 'healthy.000') %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()

3.2.3 Visualization

Open code

## plot
plotos <- data_long %>% 
  filter(!is.na(il1_b_value)) %>% 
  mutate(group = factor(group, levels = 
                          c("no_rejection", 
                            "acute_rejection",
                            "healthy")),
         time_point = factor(time_point)) %>% 
  mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>% 
  mutate(group_sep = interaction(group, time_point)) %>% 
  
## plotting
  ggplot(aes(x=time_point, y=il1_b_value)) +
  geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA, 
               size=0.5,position = position_dodge(width = 0.84), color='black') +
  labs(x = "Time post TX", y = expression("IL-1" * beta * " (pg/ml)")) +
  scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
                              "Biopsy", "Healthy")) +
  theme(axis.text = element_text(size = 11),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 12)) +
  scale_fill_manual(values = c("no_rejection" = "skyblue2",
                               "acute_rejection" = "coral",
                               "healthy" = "darkolivegreen3"), 
                    labels = c("no rejection", "acute rejection", "healthy")) +
  coord_cartesian(ylim = c(0, 28)) +
  annotate("text", x = 6, y = 10, label = "#", size = 7) 

pair_rest$y <- c(19, 17, 19, 26.5, 21, 25)

plotos <- plotos + geom_segment(data = pair_rest, 
                                aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
                                color = pair_rest$cole,
                                linetype = "solid", size = 0.4) 
plot_il1_b <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+0.5, 
                            label = pair_rest$star, size = 7, color=pair_rest$cole) 
plot_il1_b

3.3 il1_ra

3.3.1 Mixed model

Open code

## filtering data
model_data <- data_long %>% 
  filter(!is.na(il1_ra_value))

## model fit
model_il1_ra <-  lmer(log10(il1_ra_value) ~ 
                       group_sep +
                       (1|id), data = model_data)

## diagnostic plot of model
plot(model_il1_ra)

Open code

## anova result od model
Anova(model_il1_ra, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_ra_value)
##                F Df Df.res    Pr(>F)    
## group_sep 17.915  9 414.32 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## model summary
summary(model_il1_ra)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_ra_value) ~ group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 325.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.3924 -0.5919 -0.1109  0.4126  3.4509 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.02780  0.1667  
##  Residual             0.08388  0.2896  
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      3.17638    0.05311  59.805
## group_sephealthy.000            -0.57937    0.07175  -8.075
## group_sepno_rejection.000        0.11718    0.06332   1.851
## group_sepacute_rejection.007    -0.12171    0.06428  -1.893
## group_sepno_rejection.007       -0.06739    0.06332  -1.064
## group_sepacute_rejection.090    -0.11194    0.06528  -1.715
## group_sepno_rejection.090       -0.06636    0.07428  -0.893
## group_sepacute_rejection.365    -0.28796    0.06841  -4.209
## group_sepno_rejection.365       -0.17831    0.07428  -2.400
## group_sepacute_rejection.biopsy -0.17195    0.07310  -2.352
## 
## Correlation of Fixed Effects:
##              (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000  -0.740                                                      
## grp_sp_.000  -0.839  0.621                                               
## grp_spc_.007 -0.639  0.473  0.536                                        
## grp_spn_.007 -0.839  0.621  0.777   0.536                                
## grp_spc_.090 -0.632  0.468  0.530   0.520        0.530                   
## grp_spn_.090 -0.715  0.529  0.663   0.457        0.663        0.452      
## grp_spc_.365 -0.605  0.448  0.508   0.498        0.508        0.492      
## grp_spn_.365 -0.715  0.529  0.663   0.457        0.663        0.452      
## grp_spct_r.  -0.563  0.417  0.472   0.463        0.472        0.459      
##              grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000                                        
## grp_sp_.000                                        
## grp_spc_.007                                       
## grp_spn_.007                                       
## grp_spc_.090                                       
## grp_spn_.090                                       
## grp_spc_.365  0.433                                
## grp_spn_.365  0.610        0.433                   
## grp_spct_r.   0.403        0.433        0.403

3.3.2 Pairwise comparison

3.3.2.1 List of contrasts

List of contrast names and left and right coordinates on the X axis

Open code
contrast <- c(

  ## healthy vs. others
  "healthy vs. no rejection D_0",
  "healthy vs. no rejection D_7",
  "healthy vs. no rejection D_90",
  "healthy vs. no rejection D_365",
  "healthy vs. acute rejection D_0",
  "healthy vs. acute rejection D_7",
  "healthy vs. acute rejection D_90",
  "healthy vs. acute rejection D_365",
  "healthy vs. acute rejection D_biopsy",


  ## between-subjects comparisons
  "no rejection D_0 vs. acute rejection D_0",
  "no rejection D_7 vs. acute rejection D_7",
  "no rejection D_90 vs. acute rejection D_90",
  "no rejection D_365 vs. acute rejection D_365",
  "no rejection D_0 vs. acute rejection D_biopsy",
  "no rejection D_7 vs. acute rejection D_biopsy",
  "no rejection D_90 vs. acute rejection D_biopsy",
  "no rejection D_365 vs. acute rejection D_biopsy",

  ## within-subjects comparisons
  ### no rejection
  "no rejection D_0 vs. no rejection D_7",
  "no rejection D_0 vs. no rejection D_90",
  "no rejection D_0 vs. no rejection D_365",
  "no rejection D_7 vs. no rejection D_90",
  "no rejection D_7 vs. no rejection D_365",
  "no rejection D_90 vs. no rejection D_365",

    ### acute rejection
  "acute rejection D_0 vs. acute rejection D_7",
  "acute rejection D_0 vs. acute rejection D_90",
  "acute rejection D_0 vs. acute rejection D_365",
  "acute rejection D_7 vs. acute rejection D_90",
  "acute rejection D_7 vs. acute rejection D_365",
  "acute rejection D_90 vs. acute rejection D_365",
  "acute rejection D_0 vs. acute rejection D_biopsy",
  "acute rejection D_7 vs. acute rejection D_biopsy",
  "acute rejection D_90 vs. acute rejection D_biopsy",
  "acute rejection D_365 vs. acute rejection D_biopsy"
  )
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
              rep(0.8, 3), rep(1.8, 2), 2.8,
              rep(1.2, 3), rep(2.2, 2), 3.2,
              seq(1.2, 4.2, 1))

xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
              5, seq(1.2, 4.2, 1), rep(5, 4),
              seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
              2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
              rep(5, 4))

anot_coord <- (xcoord_1+xcoord_2)/2

pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
                         anot_coord)

pair_table <- pair_table %>% 
  mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>% 
  mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
  mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>% 
  mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
  
  pair_table$p <- NA

3.3.2.2 Between-subjects comparisons

Open code
set.seed(16)
for (i in 1:17){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i])|
           group_sep == paste(pair_table$term2[i]))
  
  g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
  g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
  
  pair_table$p[i] <- wilcox.test(g1$il1_ra_value, g2$il1_ra_value, paired=FALSE)$p.value
}

3.3.2.3 Within-subject comparisons

Open code
set.seed(16)
for (i in 18:33){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i]) |
           group_sep == paste(pair_table$term2[i])) %>% 
    select(id, group_sep, il1_ra_value) %>% 
    pivot_wider(names_from = group_sep, values_from = il1_ra_value) %>% 
    filter(complete.cases(.)) %>% data.frame()
  
  g1 <- dat_pair[,2]
  g2 <- dat_pair[,3]
  
  pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}

pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')

3.3.2.4 Resulting table

Open code
## the whole `pair table`
pair_table
##                                           contrast xcoord_1 xcoord_2 anot_coord
## 1                     healthy vs. no rejection D_0      6.0      0.8        3.4
## 2                     healthy vs. no rejection D_7      6.0      1.8        3.9
## 3                    healthy vs. no rejection D_90      6.0      2.8        4.4
## 4                   healthy vs. no rejection D_365      6.0      3.8        4.9
## 5                  healthy vs. acute rejection D_0      6.0      1.2        3.6
## 6                  healthy vs. acute rejection D_7      6.0      2.2        4.1
## 7                 healthy vs. acute rejection D_90      6.0      3.2        4.6
## 8                healthy vs. acute rejection D_365      6.0      4.2        5.1
## 9             healthy vs. acute rejection D_biopsy      6.0      5.0        5.5
## 10        no rejection D_0 vs. acute rejection D_0      0.8      1.2        1.0
## 11        no rejection D_7 vs. acute rejection D_7      1.8      2.2        2.0
## 12      no rejection D_90 vs. acute rejection D_90      2.8      3.2        3.0
## 13    no rejection D_365 vs. acute rejection D_365      3.8      4.2        4.0
## 16  no rejection D_90 vs. acute rejection D_biopsy      2.8      5.0        3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy      3.8      5.0        4.4
## 18           no rejection D_0 vs. no rejection D_7      0.8      1.8        1.3
## 21          no rejection D_7 vs. no rejection D_90      1.8      2.8        2.3
## 23        no rejection D_90 vs. no rejection D_365      2.8      3.8        3.3
## 24     acute rejection D_0 vs. acute rejection D_7      1.2      2.2        1.7
## 27    acute rejection D_7 vs. acute rejection D_90      2.2      3.2        2.7
## 29  acute rejection D_90 vs. acute rejection D_365      3.2      4.2        3.7
##                  term1                  term2            p  cole          fdr
## 1          healthy.000       no_rejection.000 6.278863e-18 black 1.318561e-16
## 2          healthy.000       no_rejection.007 4.243487e-16 black 4.455661e-15
## 3          healthy.000       no_rejection.090 1.424391e-11 black 9.970737e-11
## 4          healthy.000       no_rejection.365 6.321421e-09 black 1.896426e-08
## 5          healthy.000    acute_rejection.000 1.340189e-10 black 7.035992e-10
## 6          healthy.000    acute_rejection.007 2.975331e-10 black 1.249639e-09
## 7          healthy.000    acute_rejection.090 2.122261e-09 black 7.427913e-09
## 8          healthy.000    acute_rejection.365 9.383021e-06 black 2.189372e-05
## 9          healthy.000 acute_rejection.biopsy 1.355816e-07 black 3.559017e-07
## 10    no_rejection.000    acute_rejection.000 1.894461e-01 black 2.340216e-01
## 11    no_rejection.007    acute_rejection.007 6.256756e-02 black 9.581960e-02
## 12    no_rejection.090    acute_rejection.090 2.120715e-01 black 2.474167e-01
## 13    no_rejection.365    acute_rejection.365 6.387973e-02 black 9.581960e-02
## 16    no_rejection.090 acute_rejection.biopsy 1.110122e-01 black 1.554170e-01
## 17    no_rejection.365 acute_rejection.biopsy 8.746157e-01 black 8.746157e-01
## 18    no_rejection.000       no_rejection.007 6.409870e-05 blue3 1.346073e-04
## 21    no_rejection.007       no_rejection.090 5.571619e-01 blue3 6.158106e-01
## 23    no_rejection.090       no_rejection.365 2.581426e-02 blue3 4.517495e-02
## 24 acute_rejection.000    acute_rejection.007 1.422637e-01  red4 1.867211e-01
## 27 acute_rejection.007    acute_rejection.090 7.546941e-01  red4 7.924288e-01
## 29 acute_rejection.090    acute_rejection.365 2.400159e-02  red4 4.517495e-02

## only significant differences
pair_table %>% 
  filter(p<0.05) %>% 
  select(contrast, p, fdr) %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()
##                                          contrast            p          fdr
## 1                    healthy vs. no rejection D_0 6.278863e-18 1.318561e-16
## 2                    healthy vs. no rejection D_7 4.243487e-16 4.455661e-15
## 3                   healthy vs. no rejection D_90 1.424391e-11 9.970737e-11
## 4                  healthy vs. no rejection D_365 6.321421e-09 1.896426e-08
## 5                 healthy vs. acute rejection D_0 1.340189e-10 7.035992e-10
## 6                 healthy vs. acute rejection D_7 2.975331e-10 1.249639e-09
## 7                healthy vs. acute rejection D_90 2.122261e-09 7.427913e-09
## 8               healthy vs. acute rejection D_365 9.383021e-06 2.189372e-05
## 9            healthy vs. acute rejection D_biopsy 1.355816e-07 3.559017e-07
## 10          no rejection D_0 vs. no rejection D_7 6.409870e-05 1.346073e-04
## 11       no rejection D_90 vs. no rejection D_365 2.581426e-02 4.517495e-02
## 12 acute rejection D_90 vs. acute rejection D_365 2.400159e-02 4.517495e-02
##    star
## 1   ***
## 2   ***
## 3   ***
## 4   ***
## 5   ***
## 6   ***
## 7   ***
## 8   ***
## 9   ***
## 10  ***
## 11    *
## 12    *

pair_rest <- pair_table %>% 
  filter(p<0.05,
         term1 != 'healthy.000') %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()

3.3.3 Visualization

Open code

## plot
plotos <- data_long %>% 
  filter(!is.na(il1_ra_value)) %>% 
  mutate(group = factor(group, levels = 
                          c("no_rejection", 
                            "acute_rejection",
                            "healthy")),
         time_point = factor(time_point)) %>% 
  mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>% 
  mutate(group_sep = interaction(group, time_point)) %>% 
  
## plotting
  ggplot(aes(x=time_point, y=il1_ra_value)) +
  geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA, 
               size=0.5,position = position_dodge(width = 0.84), color='black') +
  labs(x = "Time post TX", y = "IL-1 RA (pg/ml)") +
  scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
                              "Biopsy", "Healthy")) +
  theme(axis.text = element_text(size = 11),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 12)) +
  scale_fill_manual(values = c("no_rejection" = "skyblue2",
                               "acute_rejection" = "coral",
                               "healthy" = "darkolivegreen3"), 
                    labels = c("no rejection", "acute rejection", "healthy")) +
  coord_cartesian(ylim = c(0, 7000)) +
  annotate("text", x = 6, y = 2000, label = "#", size = 7) 

pair_rest$y <- c(6900, 4000, 3600)

plotos <- plotos + geom_segment(data = pair_rest, 
                                aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
                                color = pair_rest$cole,
                                linetype = "solid", size = 0.4) 
plot_il1_ra <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+50, 
                            label = pair_rest$star, size = 7, color=pair_rest$cole) 
plot_il1_ra

3.4 il18

3.4.1 Mixed model

Open code

## filtering data
model_data <- data_long %>% 
  filter(!is.na(il18_value))

## model fit
model_il18 <-  lmer(log10(il18_value) ~ 
                       group_sep +
                       (1|id), data = model_data)

## diagnostic plot of model
plot(model_il18)

Open code

## anova result od model
Anova(model_il18, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_value)
##                F Df Df.res    Pr(>F)    
## group_sep 23.496  9 412.76 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## model summary
summary(model_il18)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_value) ~ group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 290.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4547 -0.5610 -0.0786  0.4362  3.2320 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.03004  0.1733  
##  Residual             0.07579  0.2753  
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                  Estimate Std. Error t value
## (Intercept)                      2.706285   0.051611  52.436
## group_sephealthy.000            -0.433108   0.069774  -6.207
## group_sepno_rejection.000        0.181152   0.061559   2.943
## group_sepacute_rejection.007    -0.250706   0.061125  -4.102
## group_sepno_rejection.007       -0.173250   0.061559  -2.814
## group_sepacute_rejection.090     0.133694   0.062099   2.153
## group_sepno_rejection.090        0.124231   0.071958   1.726
## group_sepacute_rejection.365     0.050040   0.065106   0.769
## group_sepno_rejection.365       -0.003025   0.071958  -0.042
## group_sepacute_rejection.biopsy  0.127269   0.069595   1.829
## 
## Correlation of Fixed Effects:
##              (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000  -0.740                                                      
## grp_sp_.000  -0.838  0.620                                               
## grp_spc_.007 -0.626  0.463  0.525                                        
## grp_spn_.007 -0.838  0.620  0.787   0.525                                
## grp_spc_.090 -0.619  0.458  0.519   0.521        0.519                   
## grp_spn_.090 -0.717  0.531  0.673   0.449        0.673        0.444      
## grp_spc_.365 -0.593  0.438  0.497   0.498        0.497        0.492      
## grp_spn_.365 -0.717  0.531  0.673   0.449        0.673        0.444      
## grp_spct_r.  -0.551  0.407  0.462   0.462        0.462        0.458      
##              grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000                                        
## grp_sp_.000                                        
## grp_spc_.007                                       
## grp_spn_.007                                       
## grp_spc_.090                                       
## grp_spn_.090                                       
## grp_spc_.365  0.425                                
## grp_spn_.365  0.625        0.425                   
## grp_spct_r.   0.395        0.432        0.395

3.4.2 Pairwise comparison

3.4.2.1 List of contrasts

List of contrast names and left and right coordinates on the X axis

Open code
contrast <- c(

  ## healthy vs. others
  "healthy vs. no rejection D_0",
  "healthy vs. no rejection D_7",
  "healthy vs. no rejection D_90",
  "healthy vs. no rejection D_365",
  "healthy vs. acute rejection D_0",
  "healthy vs. acute rejection D_7",
  "healthy vs. acute rejection D_90",
  "healthy vs. acute rejection D_365",
  "healthy vs. acute rejection D_biopsy",


  ## between-subjects comparisons
  "no rejection D_0 vs. acute rejection D_0",
  "no rejection D_7 vs. acute rejection D_7",
  "no rejection D_90 vs. acute rejection D_90",
  "no rejection D_365 vs. acute rejection D_365",
  "no rejection D_0 vs. acute rejection D_biopsy",
  "no rejection D_7 vs. acute rejection D_biopsy",
  "no rejection D_90 vs. acute rejection D_biopsy",
  "no rejection D_365 vs. acute rejection D_biopsy",

  ## within-subjects comparisons
  ### no rejection
  "no rejection D_0 vs. no rejection D_7",
  "no rejection D_0 vs. no rejection D_90",
  "no rejection D_0 vs. no rejection D_365",
  "no rejection D_7 vs. no rejection D_90",
  "no rejection D_7 vs. no rejection D_365",
  "no rejection D_90 vs. no rejection D_365",

    ### acute rejection
  "acute rejection D_0 vs. acute rejection D_7",
  "acute rejection D_0 vs. acute rejection D_90",
  "acute rejection D_0 vs. acute rejection D_365",
  "acute rejection D_7 vs. acute rejection D_90",
  "acute rejection D_7 vs. acute rejection D_365",
  "acute rejection D_90 vs. acute rejection D_365",
  "acute rejection D_0 vs. acute rejection D_biopsy",
  "acute rejection D_7 vs. acute rejection D_biopsy",
  "acute rejection D_90 vs. acute rejection D_biopsy",
  "acute rejection D_365 vs. acute rejection D_biopsy"
  )
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
              rep(0.8, 3), rep(1.8, 2), 2.8,
              rep(1.2, 3), rep(2.2, 2), 3.2,
              seq(1.2, 4.2, 1))

xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
              5, seq(1.2, 4.2, 1), rep(5, 4),
              seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
              2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
              rep(5, 4))

anot_coord <- (xcoord_1+xcoord_2)/2

pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
                         anot_coord)

pair_table <- pair_table %>% 
  mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>% 
  mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
  mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>% 
  mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
  
  pair_table$p <- NA

3.4.2.2 Between-subjects comparisons

Open code
set.seed(16)
for (i in 1:17){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i])|
           group_sep == paste(pair_table$term2[i]))
  
  g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
  g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
  
  pair_table$p[i] <- wilcox.test(g1$il18_value, g2$il18_value, paired=FALSE)$p.value
}

3.4.2.3 Within-subject comparisons

Open code
set.seed(16)
for (i in 18:33){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i]) |
           group_sep == paste(pair_table$term2[i])) %>% 
    select(id, group_sep, il18_value) %>% 
    pivot_wider(names_from = group_sep, values_from = il18_value) %>% 
    filter(complete.cases(.)) %>% data.frame()
  
  g1 <- dat_pair[,2]
  g2 <- dat_pair[,3]
  
  pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}

pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')

3.4.2.4 Resulting table

Open code
## the whole `pair table`
pair_table
##                                           contrast xcoord_1 xcoord_2 anot_coord
## 1                     healthy vs. no rejection D_0      6.0      0.8        3.4
## 2                     healthy vs. no rejection D_7      6.0      1.8        3.9
## 3                    healthy vs. no rejection D_90      6.0      2.8        4.4
## 4                   healthy vs. no rejection D_365      6.0      3.8        4.9
## 5                  healthy vs. acute rejection D_0      6.0      1.2        3.6
## 6                  healthy vs. acute rejection D_7      6.0      2.2        4.1
## 7                 healthy vs. acute rejection D_90      6.0      3.2        4.6
## 8                healthy vs. acute rejection D_365      6.0      4.2        5.1
## 9             healthy vs. acute rejection D_biopsy      6.0      5.0        5.5
## 10        no rejection D_0 vs. acute rejection D_0      0.8      1.2        1.0
## 11        no rejection D_7 vs. acute rejection D_7      1.8      2.2        2.0
## 12      no rejection D_90 vs. acute rejection D_90      2.8      3.2        3.0
## 13    no rejection D_365 vs. acute rejection D_365      3.8      4.2        4.0
## 16  no rejection D_90 vs. acute rejection D_biopsy      2.8      5.0        3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy      3.8      5.0        4.4
## 18           no rejection D_0 vs. no rejection D_7      0.8      1.8        1.3
## 21          no rejection D_7 vs. no rejection D_90      1.8      2.8        2.3
## 23        no rejection D_90 vs. no rejection D_365      2.8      3.8        3.3
## 24     acute rejection D_0 vs. acute rejection D_7      1.2      2.2        1.7
## 27    acute rejection D_7 vs. acute rejection D_90      2.2      3.2        2.7
## 29  acute rejection D_90 vs. acute rejection D_365      3.2      4.2        3.7
##                  term1                  term2            p  cole          fdr
## 1          healthy.000       no_rejection.000 1.671896e-16 black 3.510982e-15
## 2          healthy.000       no_rejection.007 8.598524e-08 black 2.006322e-07
## 3          healthy.000       no_rejection.090 2.667314e-13 black 1.867120e-12
## 4          healthy.000       no_rejection.365 2.761390e-10 black 8.284170e-10
## 5          healthy.000    acute_rejection.000 2.124780e-06 black 4.056397e-06
## 6          healthy.000    acute_rejection.007 4.471962e-02 black 5.524189e-02
## 7          healthy.000    acute_rejection.090 8.245243e-13 black 4.328752e-12
## 8          healthy.000    acute_rejection.365 1.675587e-12 black 7.037463e-12
## 9          healthy.000 acute_rejection.biopsy 4.551076e-11 black 1.592877e-10
## 10    no_rejection.000    acute_rejection.000 2.664446e-02 black 3.497086e-02
## 11    no_rejection.007    acute_rejection.007 1.305168e-01 black 1.522696e-01
## 12    no_rejection.090    acute_rejection.090 7.088759e-01 black 7.088759e-01
## 13    no_rejection.365    acute_rejection.365 1.722459e-02 black 2.411443e-02
## 16    no_rejection.090 acute_rejection.biopsy 3.083806e-01 black 3.237997e-01
## 17    no_rejection.365 acute_rejection.biopsy 2.268561e-03 black 3.664598e-03
## 18    no_rejection.000       no_rejection.007 1.006443e-13 blue3 1.056765e-12
## 21    no_rejection.007       no_rejection.090 1.011576e-07 blue3 2.124310e-07
## 23    no_rejection.090       no_rejection.365 6.539723e-03 blue3 9.809584e-03
## 24 acute_rejection.000    acute_rejection.007 4.081365e-05  red4 7.142390e-05
## 27 acute_rejection.007    acute_rejection.090 3.163223e-09  red4 8.303459e-09
## 29 acute_rejection.090    acute_rejection.365 2.312297e-01  red4 2.555696e-01

## only significant differences
pair_table %>% 
  filter(p<0.05) %>% 
  select(contrast, p, fdr) %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()
##                                           contrast            p          fdr
## 1                     healthy vs. no rejection D_0 1.671896e-16 3.510982e-15
## 2                     healthy vs. no rejection D_7 8.598524e-08 2.006322e-07
## 3                    healthy vs. no rejection D_90 2.667314e-13 1.867120e-12
## 4                   healthy vs. no rejection D_365 2.761390e-10 8.284170e-10
## 5                  healthy vs. acute rejection D_0 2.124780e-06 4.056397e-06
## 6                  healthy vs. acute rejection D_7 4.471962e-02 5.524189e-02
## 7                 healthy vs. acute rejection D_90 8.245243e-13 4.328752e-12
## 8                healthy vs. acute rejection D_365 1.675587e-12 7.037463e-12
## 9             healthy vs. acute rejection D_biopsy 4.551076e-11 1.592877e-10
## 10        no rejection D_0 vs. acute rejection D_0 2.664446e-02 3.497086e-02
## 11    no rejection D_365 vs. acute rejection D_365 1.722459e-02 2.411443e-02
## 12 no rejection D_365 vs. acute rejection D_biopsy 2.268561e-03 3.664598e-03
## 13           no rejection D_0 vs. no rejection D_7 1.006443e-13 1.056765e-12
## 14          no rejection D_7 vs. no rejection D_90 1.011576e-07 2.124310e-07
## 15        no rejection D_90 vs. no rejection D_365 6.539723e-03 9.809584e-03
## 16     acute rejection D_0 vs. acute rejection D_7 4.081365e-05 7.142390e-05
## 17    acute rejection D_7 vs. acute rejection D_90 3.163223e-09 8.303459e-09
##    star
## 1   ***
## 2   ***
## 3   ***
## 4   ***
## 5   ***
## 6     *
## 7   ***
## 8   ***
## 9   ***
## 10    *
## 11    *
## 12   **
## 13  ***
## 14  ***
## 15   **
## 16  ***
## 17  ***

pair_rest <- pair_table %>% 
  filter(p<0.05,
         term1 != 'healthy.000') %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()

3.4.3 Visualization

Open code

## plot
plotos <- data_long %>% 
  filter(!is.na(il18_value)) %>% 
  mutate(group = factor(group, levels = 
                          c("no_rejection", 
                            "acute_rejection",
                            "healthy")),
         time_point = factor(time_point)) %>% 
  mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>% 
  mutate(group_sep = interaction(group, time_point)) %>% 
  
## plotting
  ggplot(aes(x=time_point, y=il18_value)) +
  geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA, 
               size=0.5,position = position_dodge(width = 0.84), color='black') +
  labs(x = "Time post TX", y = "IL-18 (pg/ml)") +
  scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
                              "Biopsy", "Healthy")) +
  theme(axis.text = element_text(size = 11),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 12)) +
  scale_fill_manual(values = c("no_rejection" = "skyblue2",
                               "acute_rejection" = "coral",
                               "healthy" = "darkolivegreen3"), 
                    labels = c("no rejection", "acute rejection", "healthy")) +
  coord_cartesian(ylim = c(0, 3000)) +
  annotate("text", x = 6, y = 700, label = "#", size = 7) 

pair_rest$y <- c(2700, 1000, 1700, 2900, 2300, 2100, 1900, 1700)

plotos <- plotos + geom_segment(data = pair_rest, 
                                aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
                                color = pair_rest$cole,
                                linetype = "solid", size = 0.4) 
plot_il18 <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+50, 
                            label = pair_rest$star, size = 7, color=pair_rest$cole) 
plot_il18

3.5 il18_bp

3.5.1 Mixed model

Open code

## filtering data
model_data <- data_long %>% 
  filter(!is.na(il18_bp_value))

## model fit
model_il18_bp <-  lmer(log10(il18_bp_value) ~ 
                       group_sep +
                       (1|id), data = model_data)

## diagnostic plot of model
plot(model_il18_bp)

Open code

## anova result od model
Anova(model_il18_bp, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_bp_value)
##                F Df Df.res    Pr(>F)    
## group_sep 8.4695  9  489.8 8.293e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## model summary
summary(model_il18_bp)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_bp_value) ~ group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: -484.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.0578  -0.3070   0.0384   0.3462   4.8019 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.006046 0.07775 
##  Residual             0.018843 0.13727 
## Number of obs: 577, groups:  id, 169
## 
## Fixed effects:
##                                  Estimate Std. Error t value
## (Intercept)                      3.614819   0.024804 145.735
## group_sephealthy.000            -0.165875   0.037658  -4.405
## group_sepno_rejection.000        0.004244   0.029707   0.143
## group_sepacute_rejection.007    -0.061525   0.030101  -2.044
## group_sepno_rejection.007       -0.046670   0.029665  -1.573
## group_sepacute_rejection.090    -0.121220   0.030711  -3.947
## group_sepno_rejection.090       -0.109015   0.029982  -3.636
## group_sepacute_rejection.365    -0.092007   0.031914  -2.883
## group_sepno_rejection.365       -0.121527   0.030184  -4.026
## group_sepacute_rejection.biopsy -0.053131   0.034328  -1.548
## 
## Correlation of Fixed Effects:
##              (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000  -0.659                                                      
## grp_sp_.000  -0.835  0.550                                               
## grp_spc_.007 -0.640  0.422  0.534                                        
## grp_spn_.007 -0.836  0.551  0.771   0.535                                
## grp_spc_.090 -0.628  0.414  0.524   0.517        0.525                   
## grp_spn_.090 -0.827  0.545  0.763   0.529        0.764        0.520      
## grp_spc_.365 -0.606  0.399  0.506   0.499        0.507        0.489      
## grp_spn_.365 -0.822  0.541  0.758   0.526        0.759        0.516      
## grp_spct_r.  -0.556  0.366  0.464   0.458        0.465        0.449      
##              grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000                                        
## grp_sp_.000                                        
## grp_spc_.007                                       
## grp_spn_.007                                       
## grp_spc_.090                                       
## grp_spn_.090                                       
## grp_spc_.365  0.501                                
## grp_spn_.365  0.752        0.498                   
## grp_spct_r.   0.460        0.428        0.457

3.5.2 Pairwise comparison

3.5.2.1 List of contrasts

List of contrast names and left and right coordinates on the X axis

Open code
contrast <- c(

  ## healthy vs. others
  "healthy vs. no rejection D_0",
  "healthy vs. no rejection D_7",
  "healthy vs. no rejection D_90",
  "healthy vs. no rejection D_365",
  "healthy vs. acute rejection D_0",
  "healthy vs. acute rejection D_7",
  "healthy vs. acute rejection D_90",
  "healthy vs. acute rejection D_365",
  "healthy vs. acute rejection D_biopsy",


  ## between-subjects comparisons
  "no rejection D_0 vs. acute rejection D_0",
  "no rejection D_7 vs. acute rejection D_7",
  "no rejection D_90 vs. acute rejection D_90",
  "no rejection D_365 vs. acute rejection D_365",
  "no rejection D_0 vs. acute rejection D_biopsy",
  "no rejection D_7 vs. acute rejection D_biopsy",
  "no rejection D_90 vs. acute rejection D_biopsy",
  "no rejection D_365 vs. acute rejection D_biopsy",

  ## within-subjects comparisons
  ### no rejection
  "no rejection D_0 vs. no rejection D_7",
  "no rejection D_0 vs. no rejection D_90",
  "no rejection D_0 vs. no rejection D_365",
  "no rejection D_7 vs. no rejection D_90",
  "no rejection D_7 vs. no rejection D_365",
  "no rejection D_90 vs. no rejection D_365",

    ### acute rejection
  "acute rejection D_0 vs. acute rejection D_7",
  "acute rejection D_0 vs. acute rejection D_90",
  "acute rejection D_0 vs. acute rejection D_365",
  "acute rejection D_7 vs. acute rejection D_90",
  "acute rejection D_7 vs. acute rejection D_365",
  "acute rejection D_90 vs. acute rejection D_365",
  "acute rejection D_0 vs. acute rejection D_biopsy",
  "acute rejection D_7 vs. acute rejection D_biopsy",
  "acute rejection D_90 vs. acute rejection D_biopsy",
  "acute rejection D_365 vs. acute rejection D_biopsy"
  )
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
              rep(0.8, 3), rep(1.8, 2), 2.8,
              rep(1.2, 3), rep(2.2, 2), 3.2,
              seq(1.2, 4.2, 1))

xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
              5, seq(1.2, 4.2, 1), rep(5, 4),
              seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
              2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
              rep(5, 4))

anot_coord <- (xcoord_1+xcoord_2)/2

pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
                         anot_coord)

pair_table <- pair_table %>% 
  mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>% 
  mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
  mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>% 
  mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
  
  pair_table$p <- NA

3.5.2.2 Between-subjects comparisons

Open code
set.seed(16)
for (i in 1:17){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i])|
           group_sep == paste(pair_table$term2[i]))
  
  g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
  g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
  
  pair_table$p[i] <- wilcox.test(g1$il18_bp_value, g2$il18_bp_value, paired=FALSE)$p.value
}

3.5.2.3 Within-subject comparisons

Open code
set.seed(16)
for (i in 18:33){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i]) |
           group_sep == paste(pair_table$term2[i])) %>% 
    select(id, group_sep, il18_bp_value) %>% 
    pivot_wider(names_from = group_sep, values_from = il18_bp_value) %>% 
    filter(complete.cases(.)) %>% data.frame()
  
  g1 <- dat_pair[,2]
  g2 <- dat_pair[,3]
  
  pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}

pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')

3.5.2.4 Resulting table

Open code
## the whole `pair table`
pair_table
##                                           contrast xcoord_1 xcoord_2 anot_coord
## 1                     healthy vs. no rejection D_0      6.0      0.8        3.4
## 2                     healthy vs. no rejection D_7      6.0      1.8        3.9
## 3                    healthy vs. no rejection D_90      6.0      2.8        4.4
## 4                   healthy vs. no rejection D_365      6.0      3.8        4.9
## 5                  healthy vs. acute rejection D_0      6.0      1.2        3.6
## 6                  healthy vs. acute rejection D_7      6.0      2.2        4.1
## 7                 healthy vs. acute rejection D_90      6.0      3.2        4.6
## 8                healthy vs. acute rejection D_365      6.0      4.2        5.1
## 9             healthy vs. acute rejection D_biopsy      6.0      5.0        5.5
## 10        no rejection D_0 vs. acute rejection D_0      0.8      1.2        1.0
## 11        no rejection D_7 vs. acute rejection D_7      1.8      2.2        2.0
## 12      no rejection D_90 vs. acute rejection D_90      2.8      3.2        3.0
## 13    no rejection D_365 vs. acute rejection D_365      3.8      4.2        4.0
## 16  no rejection D_90 vs. acute rejection D_biopsy      2.8      5.0        3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy      3.8      5.0        4.4
## 18           no rejection D_0 vs. no rejection D_7      0.8      1.8        1.3
## 21          no rejection D_7 vs. no rejection D_90      1.8      2.8        2.3
## 23        no rejection D_90 vs. no rejection D_365      2.8      3.8        3.3
## 24     acute rejection D_0 vs. acute rejection D_7      1.2      2.2        1.7
## 27    acute rejection D_7 vs. acute rejection D_90      2.2      3.2        2.7
## 29  acute rejection D_90 vs. acute rejection D_365      3.2      4.2        3.7
##                  term1                  term2            p  cole          fdr
## 1          healthy.000       no_rejection.000 4.111120e-12 black 8.633353e-11
## 2          healthy.000       no_rejection.007 1.240228e-08 black 8.681598e-08
## 3          healthy.000       no_rejection.090 2.278600e-03 black 4.350054e-03
## 4          healthy.000       no_rejection.365 1.381490e-04 black 3.626412e-04
## 5          healthy.000    acute_rejection.000 1.711485e-09 black 1.797059e-08
## 6          healthy.000    acute_rejection.007 9.191676e-05 black 2.757503e-04
## 7          healthy.000    acute_rejection.090 4.719442e-02 black 6.607219e-02
## 8          healthy.000    acute_rejection.365 8.426347e-03 black 1.361179e-02
## 9          healthy.000 acute_rejection.biopsy 1.333362e-06 black 5.600122e-06
## 10    no_rejection.000    acute_rejection.000 3.168633e-01 black 3.696738e-01
## 11    no_rejection.007    acute_rejection.007 1.275900e-01 black 1.674619e-01
## 12    no_rejection.090    acute_rejection.090 4.543678e-01 black 5.021959e-01
## 13    no_rejection.365    acute_rejection.365 7.350379e-01 black 7.350379e-01
## 16    no_rejection.090 acute_rejection.biopsy 4.632816e-03 black 8.107427e-03
## 17    no_rejection.365 acute_rejection.biopsy 3.163257e-02 black 4.744886e-02
## 18    no_rejection.000       no_rejection.007 2.090599e-05 blue3 7.317096e-05
## 21    no_rejection.007       no_rejection.090 4.544604e-08 blue3 2.385917e-07
## 23    no_rejection.090       no_rejection.365 5.092835e-01 blue3 5.347477e-01
## 24 acute_rejection.000    acute_rejection.007 1.367483e-03  red4 2.943894e-03
## 27 acute_rejection.007    acute_rejection.090 1.401854e-03  red4 2.943894e-03
## 29 acute_rejection.090    acute_rejection.365 1.547263e-01  red4 1.911325e-01

## only significant differences
pair_table %>% 
  filter(p<0.05) %>% 
  select(contrast, p, fdr) %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()
##                                           contrast            p          fdr
## 1                     healthy vs. no rejection D_0 4.111120e-12 8.633353e-11
## 2                     healthy vs. no rejection D_7 1.240228e-08 8.681598e-08
## 3                    healthy vs. no rejection D_90 2.278600e-03 4.350054e-03
## 4                   healthy vs. no rejection D_365 1.381490e-04 3.626412e-04
## 5                  healthy vs. acute rejection D_0 1.711485e-09 1.797059e-08
## 6                  healthy vs. acute rejection D_7 9.191676e-05 2.757503e-04
## 7                 healthy vs. acute rejection D_90 4.719442e-02 6.607219e-02
## 8                healthy vs. acute rejection D_365 8.426347e-03 1.361179e-02
## 9             healthy vs. acute rejection D_biopsy 1.333362e-06 5.600122e-06
## 10  no rejection D_90 vs. acute rejection D_biopsy 4.632816e-03 8.107427e-03
## 11 no rejection D_365 vs. acute rejection D_biopsy 3.163257e-02 4.744886e-02
## 12           no rejection D_0 vs. no rejection D_7 2.090599e-05 7.317096e-05
## 13          no rejection D_7 vs. no rejection D_90 4.544604e-08 2.385917e-07
## 14     acute rejection D_0 vs. acute rejection D_7 1.367483e-03 2.943894e-03
## 15    acute rejection D_7 vs. acute rejection D_90 1.401854e-03 2.943894e-03
##    star
## 1   ***
## 2   ***
## 3    **
## 4   ***
## 5   ***
## 6   ***
## 7     *
## 8    **
## 9   ***
## 10   **
## 11    *
## 12  ***
## 13  ***
## 14   **
## 15   **

pair_rest <- pair_table %>% 
  filter(p<0.05,
         term1 != 'healthy.000') %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()

! Healthy group does not differ from others with p<0.001, but some only p<0.05

3.5.3 Visualization

Open code

## plot
plotos <- data_long %>% 
  filter(!is.na(il18_bp_value)) %>% 
  mutate(group = factor(group, levels = 
                          c("no_rejection", 
                            "acute_rejection",
                            "healthy")),
         time_point = factor(time_point)) %>% 
  mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>% 
  mutate(group_sep = interaction(group, time_point)) %>% 
  
## plotting
  ggplot(aes(x=time_point, y=il18_bp_value)) +
  geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA, 
               size=0.5,position = position_dodge(width = 0.84), color='black') +
  labs(x = "Time post TX", y = "IL-18 BP (pg/ml)") +
  scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
                              "Biopsy", "Healthy")) +
  theme(axis.text = element_text(size = 11),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 12)) +
  scale_fill_manual(values = c("no_rejection" = "skyblue2",
                               "acute_rejection" = "coral",
                               "healthy" = "darkolivegreen3"), 
                    labels = c("no rejection", "acute rejection", "healthy")) +
  coord_cartesian(ylim = c(0, 7300)) +
  annotate("text", x = 6, y = 5000, label = "#", size = 7) 

pair_rest$y <- c(5800, 5500, 7200, 6900, 6400, 6100)

plotos <- plotos + geom_segment(data = pair_rest, 
                                aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
                                color = pair_rest$cole,
                                linetype = "solid", size = 0.4) 
plot_il18_bp <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+50, 
                            label = pair_rest$star, size = 7, color=pair_rest$cole) 
plot_il18_bp

3.6 il18_free

3.6.1 Mixed model

Open code

## filtering data
model_data <- data_long %>% 
  filter(!is.na(il18_free))

## model fit
model_il18_free <-  lmer(log10(il18_free) ~ 
                       group_sep +
                       (1|id), data = model_data)

## diagnostic plot of model
plot(model_il18_free)

Open code

## anova result od model
Anova(model_il18_free, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_free)
##                F Df Df.res    Pr(>F)    
## group_sep 19.077  9 401.57 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## model summary
summary(model_il18_free)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_free) ~ group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 247.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4008 -0.5625 -0.0955  0.4128  3.3292 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.02494  0.1579  
##  Residual             0.07261  0.2695  
## Number of obs: 481, groups:  id, 169
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      2.51490    0.05020  50.101
## group_sephealthy.000            -0.33233    0.07528  -4.415
## group_sepno_rejection.000        0.16742    0.05973   2.803
## group_sepacute_rejection.007    -0.23114    0.06029  -3.834
## group_sepno_rejection.007       -0.16130    0.05964  -2.704
## group_sepacute_rejection.090     0.17291    0.06123   2.824
## group_sepno_rejection.090        0.14789    0.06978   2.119
## group_sepacute_rejection.365     0.07912    0.06415   1.233
## group_sepno_rejection.365        0.01916    0.06978   0.275
## group_sepacute_rejection.biopsy  0.14763    0.06836   2.160
## 
## Correlation of Fixed Effects:
##              (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000  -0.667                                                      
## grp_sp_.000  -0.840  0.560                                               
## grp_spc_.007 -0.643  0.429  0.540                                        
## grp_spn_.007 -0.842  0.561  0.782   0.541                                
## grp_spc_.090 -0.636  0.424  0.535   0.528        0.536                   
## grp_spn_.090 -0.719  0.480  0.669   0.463        0.669        0.458      
## grp_spc_.365 -0.610  0.407  0.513   0.506        0.514        0.500      
## grp_spn_.365 -0.719  0.480  0.669   0.463        0.669        0.458      
## grp_spct_r.  -0.566  0.378  0.476   0.469        0.477        0.465      
##              grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000                                        
## grp_sp_.000                                        
## grp_spc_.007                                       
## grp_spn_.007                                       
## grp_spc_.090                                       
## grp_spn_.090                                       
## grp_spc_.365  0.439                                
## grp_spn_.365  0.618        0.439                   
## grp_spct_r.   0.407        0.439        0.407

3.6.2 Pairwise comparison

3.6.2.1 List of contrasts

List of contrast names and left and right coordinates on the X axis

Open code
contrast <- c(

  ## healthy vs. others
  "healthy vs. no rejection D_0",
  "healthy vs. no rejection D_7",
  "healthy vs. no rejection D_90",
  "healthy vs. no rejection D_365",
  "healthy vs. acute rejection D_0",
  "healthy vs. acute rejection D_7",
  "healthy vs. acute rejection D_90",
  "healthy vs. acute rejection D_365",
  "healthy vs. acute rejection D_biopsy",


  ## between-subjects comparisons
  "no rejection D_0 vs. acute rejection D_0",
  "no rejection D_7 vs. acute rejection D_7",
  "no rejection D_90 vs. acute rejection D_90",
  "no rejection D_365 vs. acute rejection D_365",
  "no rejection D_0 vs. acute rejection D_biopsy",
  "no rejection D_7 vs. acute rejection D_biopsy",
  "no rejection D_90 vs. acute rejection D_biopsy",
  "no rejection D_365 vs. acute rejection D_biopsy",

  ## within-subjects comparisons
  ### no rejection
  "no rejection D_0 vs. no rejection D_7",
  "no rejection D_0 vs. no rejection D_90",
  "no rejection D_0 vs. no rejection D_365",
  "no rejection D_7 vs. no rejection D_90",
  "no rejection D_7 vs. no rejection D_365",
  "no rejection D_90 vs. no rejection D_365",

    ### acute rejection
  "acute rejection D_0 vs. acute rejection D_7",
  "acute rejection D_0 vs. acute rejection D_90",
  "acute rejection D_0 vs. acute rejection D_365",
  "acute rejection D_7 vs. acute rejection D_90",
  "acute rejection D_7 vs. acute rejection D_365",
  "acute rejection D_90 vs. acute rejection D_365",
  "acute rejection D_0 vs. acute rejection D_biopsy",
  "acute rejection D_7 vs. acute rejection D_biopsy",
  "acute rejection D_90 vs. acute rejection D_biopsy",
  "acute rejection D_365 vs. acute rejection D_biopsy"
  )
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
              rep(0.8, 3), rep(1.8, 2), 2.8,
              rep(1.2, 3), rep(2.2, 2), 3.2,
              seq(1.2, 4.2, 1))

xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
              5, seq(1.2, 4.2, 1), rep(5, 4),
              seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
              2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
              rep(5, 4))

anot_coord <- (xcoord_1+xcoord_2)/2

pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
                         anot_coord)

pair_table <- pair_table %>% 
  mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>% 
  mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
  mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>% 
  mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
  
  pair_table$p <- NA

3.6.2.2 Between-subjects comparisons

Open code
set.seed(16)
for (i in 1:17){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i])|
           group_sep == paste(pair_table$term2[i]))
  
  g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
  g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
  
  pair_table$p[i] <- wilcox.test(g1$il18_free, g2$il18_free, paired=FALSE)$p.value
}

3.6.2.3 Within-subject comparisons

Open code
set.seed(16)
for (i in 18:33){
  dat_pair <- model_data %>%
    filter(group_sep == paste(pair_table$term1[i]) |
           group_sep == paste(pair_table$term2[i])) %>%
    select(id, group_sep, il18_free) %>%
    pivot_wider(names_from = group_sep, values_from = il18_free) %>%
    filter(complete.cases(.)) %>% data.frame()

  g1 <- dat_pair[,2]
  g2 <- dat_pair[,3]

  pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}

pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')

3.6.2.4 Resulting table

Open code
## the whole `pair table`
pair_table
##                                           contrast xcoord_1 xcoord_2 anot_coord
## 1                     healthy vs. no rejection D_0      6.0      0.8        3.4
## 2                     healthy vs. no rejection D_7      6.0      1.8        3.9
## 3                    healthy vs. no rejection D_90      6.0      2.8        4.4
## 4                   healthy vs. no rejection D_365      6.0      3.8        4.9
## 5                  healthy vs. acute rejection D_0      6.0      1.2        3.6
## 6                  healthy vs. acute rejection D_7      6.0      2.2        4.1
## 7                 healthy vs. acute rejection D_90      6.0      3.2        4.6
## 8                healthy vs. acute rejection D_365      6.0      4.2        5.1
## 9             healthy vs. acute rejection D_biopsy      6.0      5.0        5.5
## 10        no rejection D_0 vs. acute rejection D_0      0.8      1.2        1.0
## 11        no rejection D_7 vs. acute rejection D_7      1.8      2.2        2.0
## 12      no rejection D_90 vs. acute rejection D_90      2.8      3.2        3.0
## 13    no rejection D_365 vs. acute rejection D_365      3.8      4.2        4.0
## 16  no rejection D_90 vs. acute rejection D_biopsy      2.8      5.0        3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy      3.8      5.0        4.4
## 18           no rejection D_0 vs. no rejection D_7      0.8      1.8        1.3
## 21          no rejection D_7 vs. no rejection D_90      1.8      2.8        2.3
## 23        no rejection D_90 vs. no rejection D_365      2.8      3.8        3.3
## 24     acute rejection D_0 vs. acute rejection D_7      1.2      2.2        1.7
## 27    acute rejection D_7 vs. acute rejection D_90      2.2      3.2        2.7
## 29  acute rejection D_90 vs. acute rejection D_365      3.2      4.2        3.7
##                  term1                  term2            p  cole          fdr
## 1          healthy.000       no_rejection.000 1.195713e-10 black 4.184996e-10
## 2          healthy.000       no_rejection.007 1.125705e-04 black 2.363980e-04
## 3          healthy.000       no_rejection.090 1.142557e-13 black 2.382870e-12
## 4          healthy.000       no_rejection.365 2.897492e-09 black 7.605918e-09
## 5          healthy.000    acute_rejection.000 1.390081e-04 black 2.646890e-04
## 6          healthy.000    acute_rejection.007 2.848719e-01 black 3.148584e-01
## 7          healthy.000    acute_rejection.090 2.269400e-13 black 2.382870e-12
## 8          healthy.000    acute_rejection.365 1.488787e-12 black 1.042151e-11
## 9          healthy.000 acute_rejection.biopsy 1.051022e-10 black 4.184996e-10
## 10    no_rejection.000    acute_rejection.000 3.508526e-02 black 4.604940e-02
## 11    no_rejection.007    acute_rejection.007 9.916951e-02 black 1.225035e-01
## 12    no_rejection.090    acute_rejection.090 4.724356e-01 black 4.724356e-01
## 13    no_rejection.365    acute_rejection.365 1.561394e-02 black 2.185951e-02
## 16    no_rejection.090 acute_rejection.biopsy 3.794611e-01 black 3.984342e-01
## 17    no_rejection.365 acute_rejection.biopsy 3.232913e-03 black 5.222398e-03
## 18    no_rejection.000       no_rejection.007 6.418077e-12 blue3 3.369490e-11
## 21    no_rejection.007       no_rejection.090 3.163586e-08 blue3 7.381702e-08
## 23    no_rejection.090       no_rejection.365 5.187312e-03 blue3 7.780968e-03
## 24 acute_rejection.000    acute_rejection.007 1.512508e-04  red4 2.646890e-04
## 27 acute_rejection.007    acute_rejection.090 1.644366e-09  red4 4.933099e-09
## 29 acute_rejection.090    acute_rejection.365 1.776244e-01  red4 2.072284e-01

## only significant differences
pair_table %>% 
  filter(p<0.05) %>% 
  select(contrast, p, fdr) %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()
##                                           contrast            p          fdr
## 1                     healthy vs. no rejection D_0 1.195713e-10 4.184996e-10
## 2                     healthy vs. no rejection D_7 1.125705e-04 2.363980e-04
## 3                    healthy vs. no rejection D_90 1.142557e-13 2.382870e-12
## 4                   healthy vs. no rejection D_365 2.897492e-09 7.605918e-09
## 5                  healthy vs. acute rejection D_0 1.390081e-04 2.646890e-04
## 6                 healthy vs. acute rejection D_90 2.269400e-13 2.382870e-12
## 7                healthy vs. acute rejection D_365 1.488787e-12 1.042151e-11
## 8             healthy vs. acute rejection D_biopsy 1.051022e-10 4.184996e-10
## 9         no rejection D_0 vs. acute rejection D_0 3.508526e-02 4.604940e-02
## 10    no rejection D_365 vs. acute rejection D_365 1.561394e-02 2.185951e-02
## 11 no rejection D_365 vs. acute rejection D_biopsy 3.232913e-03 5.222398e-03
## 12           no rejection D_0 vs. no rejection D_7 6.418077e-12 3.369490e-11
## 13          no rejection D_7 vs. no rejection D_90 3.163586e-08 7.381702e-08
## 14        no rejection D_90 vs. no rejection D_365 5.187312e-03 7.780968e-03
## 15     acute rejection D_0 vs. acute rejection D_7 1.512508e-04 2.646890e-04
## 16    acute rejection D_7 vs. acute rejection D_90 1.644366e-09 4.933099e-09
##    star
## 1   ***
## 2   ***
## 3   ***
## 4   ***
## 5   ***
## 6   ***
## 7   ***
## 8   ***
## 9     *
## 10    *
## 11   **
## 12  ***
## 13  ***
## 14   **
## 15  ***
## 16  ***

pair_rest <- pair_table %>% 
  filter(p<0.05,
         term1 != 'healthy.000') %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()

3.6.3 Visualization

Open code

## plot
plotos <- data_long %>% 
  filter(!is.na(il18_free)) %>% 
  mutate(group = factor(group, levels = 
                          c("no_rejection", 
                            "acute_rejection",
                            "healthy")),
         time_point = factor(time_point)) %>% 
  mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>% 
  mutate(group_sep = interaction(group, time_point)) %>% 
  
## plotting
  ggplot(aes(x=time_point, y=il18_free)) +
  geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA, 
               size=0.5,position = position_dodge(width = 0.84), color='black') +
  labs(x = "Time post TX", y = "Free IL-18 (pg/ml)") +
  scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
                              "Biopsy", "Healthy")) +
  theme(axis.text = element_text(size = 11),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 12)) +
  scale_fill_manual(values = c("no_rejection" = "skyblue2",
                               "acute_rejection" = "coral",
                               "healthy" = "darkolivegreen3"), 
                    labels = c("no rejection", "acute rejection", "healthy")) +
  coord_cartesian(ylim = c(0, 1700)) +
  annotate("text", x = 6, y = 500, label = "#", size = 7) 

pair_rest$y <- c(1650,  700, 1060, 1720, 
                 830, 1200, 1350, 1280)

plotos <- plotos + geom_segment(data = pair_rest, 
                                aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
                                color = pair_rest$cole,
                                linetype = "solid", size = 0.4) 
plot_il18_free <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+5, 
                            label = pair_rest$star, size = 7, color=pair_rest$cole) 
plot_il18_free

3.7 il36_b

3.7.1 Mixed model

Open code

## filtering data
model_data <- data_long %>% 
  filter(!is.na(il36_b_value))

## model fit
model_il36_b <-  lmer(log10(il36_b_value) ~ 
                       group_sep +
                       (1|id), data = model_data)

## diagnostic plot of model
plot(model_il36_b)

Open code

## anova result od model
Anova(model_il36_b, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il36_b_value)
##                F Df Df.res    Pr(>F)    
## group_sep 13.599  9 411.76 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## model summary
summary(model_il36_b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il36_b_value) ~ group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 566.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3627 -0.3397  0.0080  0.4167  3.2126 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.05778  0.2404  
##  Residual             0.13031  0.3610  
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      0.70905    0.06872  10.318
## group_sephealthy.000            -0.52054    0.09296  -5.600
## group_sepno_rejection.000        0.05180    0.08200   0.632
## group_sepacute_rejection.007    -0.16912    0.08016  -2.110
## group_sepno_rejection.007       -0.05161    0.08200  -0.629
## group_sepacute_rejection.090    -0.36787    0.08146  -4.516
## group_sepno_rejection.090       -0.13766    0.09559  -1.440
## group_sepacute_rejection.365    -0.52969    0.08543  -6.200
## group_sepno_rejection.365       -0.19032    0.09559  -1.991
## group_sepacute_rejection.biopsy -0.50731    0.09134  -5.554
## 
## Correlation of Fixed Effects:
##              (Intr) g_.000 g__.000 grp_spc_.007 grp_spn_.007 grp_spc_.090
## grp_sph.000  -0.739                                                      
## grp_sp_.000  -0.838  0.620                                               
## grp_spc_.007 -0.616  0.456  0.517                                        
## grp_spn_.007 -0.838  0.620  0.794   0.517                                
## grp_spc_.090 -0.610  0.451  0.511   0.521        0.511                   
## grp_spn_.090 -0.719  0.531  0.681   0.443        0.681        0.439      
## grp_spc_.365 -0.584  0.432  0.489   0.498        0.489        0.492      
## grp_spn_.365 -0.719  0.531  0.681   0.443        0.681        0.439      
## grp_spct_r.  -0.542  0.401  0.455   0.462        0.455        0.458      
##              grp_spn_.090 grp_spc_.365 grp_spn_.365
## grp_sph.000                                        
## grp_sp_.000                                        
## grp_spc_.007                                       
## grp_spn_.007                                       
## grp_spc_.090                                       
## grp_spn_.090                                       
## grp_spc_.365  0.420                                
## grp_spn_.365  0.634        0.420                   
## grp_spct_r.   0.390        0.432        0.390

3.7.2 Pairwise comparison

3.7.2.1 List of contrasts

List of contrast names and left and right coordinates on the X axis

Open code
contrast <- c(

  ## healthy vs. others
  "healthy vs. no rejection D_0",
  "healthy vs. no rejection D_7",
  "healthy vs. no rejection D_90",
  "healthy vs. no rejection D_365",
  "healthy vs. acute rejection D_0",
  "healthy vs. acute rejection D_7",
  "healthy vs. acute rejection D_90",
  "healthy vs. acute rejection D_365",
  "healthy vs. acute rejection D_biopsy",


  ## between-subjects comparisons
  "no rejection D_0 vs. acute rejection D_0",
  "no rejection D_7 vs. acute rejection D_7",
  "no rejection D_90 vs. acute rejection D_90",
  "no rejection D_365 vs. acute rejection D_365",
  "no rejection D_0 vs. acute rejection D_biopsy",
  "no rejection D_7 vs. acute rejection D_biopsy",
  "no rejection D_90 vs. acute rejection D_biopsy",
  "no rejection D_365 vs. acute rejection D_biopsy",

  ## within-subjects comparisons
  ### no rejection
  "no rejection D_0 vs. no rejection D_7",
  "no rejection D_0 vs. no rejection D_90",
  "no rejection D_0 vs. no rejection D_365",
  "no rejection D_7 vs. no rejection D_90",
  "no rejection D_7 vs. no rejection D_365",
  "no rejection D_90 vs. no rejection D_365",

    ### acute rejection
  "acute rejection D_0 vs. acute rejection D_7",
  "acute rejection D_0 vs. acute rejection D_90",
  "acute rejection D_0 vs. acute rejection D_365",
  "acute rejection D_7 vs. acute rejection D_90",
  "acute rejection D_7 vs. acute rejection D_365",
  "acute rejection D_90 vs. acute rejection D_365",
  "acute rejection D_0 vs. acute rejection D_biopsy",
  "acute rejection D_7 vs. acute rejection D_biopsy",
  "acute rejection D_90 vs. acute rejection D_biopsy",
  "acute rejection D_365 vs. acute rejection D_biopsy"
  )
xcoord_1 <- c(rep(6, 9), seq(0.8, 3.8, 1), seq(0.8, 3.8, 1),
              rep(0.8, 3), rep(1.8, 2), 2.8,
              rep(1.2, 3), rep(2.2, 2), 3.2,
              seq(1.2, 4.2, 1))

xcoord_2 <- c(seq(0.8, 3.8, 1), seq(1.2, 4.2, 1),
              5, seq(1.2, 4.2, 1), rep(5, 4),
              seq(1.8, 3.8, 1), 2.8, 3.8, 3.8,
              2.2, 3.2, 4.2, 3.2, 4.2, 4.2,
              rep(5, 4))

anot_coord <- (xcoord_1+xcoord_2)/2

pair_table <- data.frame(contrast, xcoord_1, xcoord_2,
                         anot_coord)

pair_table <- pair_table %>% 
  mutate(term1 = if_else(xcoord_1 == 6, 'healthy.000', 'other')) %>% 
  mutate(term1 = if_else(xcoord_1 == 0.8, 'no_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.8, 'no_rejection.007', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.8, 'no_rejection.090', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 3.8, 'no_rejection.365', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 5, 'acute_rejection.biopsy', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 1.2, 'acute_rejection.000', term1)) %>% 
  mutate(term1 = if_else(xcoord_1 == 2.2, 'acute_rejection.007', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 3.2, 'acute_rejection.090', term1)) %>%
  mutate(term1 = if_else(xcoord_1 == 4.2, 'acute_rejection.365', term1)) %>%
  mutate(term2 = if_else(xcoord_2 == 6, 'healthy.000', 'other')) %>% 
  mutate(term2 = if_else(xcoord_2 == 0.8, 'no_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.8, 'no_rejection.007', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.8, 'no_rejection.090', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 3.8, 'no_rejection.365', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 5, 'acute_rejection.biopsy', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 1.2, 'acute_rejection.000', term2)) %>% 
  mutate(term2 = if_else(xcoord_2 == 2.2, 'acute_rejection.007', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 3.2, 'acute_rejection.090', term2)) %>%
  mutate(term2 = if_else(xcoord_2 == 4.2, 'acute_rejection.365', term2))
  
  pair_table$p <- NA

3.7.2.2 Between-subjects comparisons

Open code
set.seed(16)
for (i in 1:17){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i])|
           group_sep == paste(pair_table$term2[i]))
  
  g1 <- dat_pair %>% filter(group_sep == paste(pair_table$term1[i]))
  g2 <- dat_pair %>% filter(group_sep == paste(pair_table$term2[i]))
  
  pair_table$p[i] <- wilcox.test(g1$il36_b_value, g2$il36_b_value, paired=FALSE)$p.value
}

3.7.2.3 Within-subject comparisons

Open code
set.seed(16)
for (i in 18:33){
  dat_pair <- model_data %>% 
    filter(group_sep == paste(pair_table$term1[i]) |
           group_sep == paste(pair_table$term2[i])) %>% 
    select(id, group_sep, il36_b_value) %>% 
    pivot_wider(names_from = group_sep, values_from = il36_b_value) %>% 
    filter(complete.cases(.)) %>% data.frame()
  
  g1 <- dat_pair[,2]
  g2 <- dat_pair[,3]
  
  pair_table$p[i] <- wilcox.test(c(g1), c(g2), paired=TRUE)$p.value
}

pair_table <- pair_table[-c(14:15, 19:20, 22, 25:26, 28, 30:33),]
pair_table$cole <- c(rep("black", 15), rep("blue3", 3), rep("red4", 3))
pair_table$fdr <- p.adjust(pair_table$p, method = 'BH')

3.7.2.4 Resulting table

Open code
## the whole `pair table`
pair_table
##                                           contrast xcoord_1 xcoord_2 anot_coord
## 1                     healthy vs. no rejection D_0      6.0      0.8        3.4
## 2                     healthy vs. no rejection D_7      6.0      1.8        3.9
## 3                    healthy vs. no rejection D_90      6.0      2.8        4.4
## 4                   healthy vs. no rejection D_365      6.0      3.8        4.9
## 5                  healthy vs. acute rejection D_0      6.0      1.2        3.6
## 6                  healthy vs. acute rejection D_7      6.0      2.2        4.1
## 7                 healthy vs. acute rejection D_90      6.0      3.2        4.6
## 8                healthy vs. acute rejection D_365      6.0      4.2        5.1
## 9             healthy vs. acute rejection D_biopsy      6.0      5.0        5.5
## 10        no rejection D_0 vs. acute rejection D_0      0.8      1.2        1.0
## 11        no rejection D_7 vs. acute rejection D_7      1.8      2.2        2.0
## 12      no rejection D_90 vs. acute rejection D_90      2.8      3.2        3.0
## 13    no rejection D_365 vs. acute rejection D_365      3.8      4.2        4.0
## 16  no rejection D_90 vs. acute rejection D_biopsy      2.8      5.0        3.9
## 17 no rejection D_365 vs. acute rejection D_biopsy      3.8      5.0        4.4
## 18           no rejection D_0 vs. no rejection D_7      0.8      1.8        1.3
## 21          no rejection D_7 vs. no rejection D_90      1.8      2.8        2.3
## 23        no rejection D_90 vs. no rejection D_365      2.8      3.8        3.3
## 24     acute rejection D_0 vs. acute rejection D_7      1.2      2.2        1.7
## 27    acute rejection D_7 vs. acute rejection D_90      2.2      3.2        2.7
## 29  acute rejection D_90 vs. acute rejection D_365      3.2      4.2        3.7
##                  term1                  term2            p  cole          fdr
## 1          healthy.000       no_rejection.000 1.507542e-15 black 3.165839e-14
## 2          healthy.000       no_rejection.007 1.250083e-13 black 1.312587e-12
## 3          healthy.000       no_rejection.090 5.134600e-09 black 2.695665e-08
## 4          healthy.000       no_rejection.365 7.018490e-08 black 2.947766e-07
## 5          healthy.000    acute_rejection.000 1.419560e-11 black 9.936919e-11
## 6          healthy.000    acute_rejection.007 1.164937e-06 black 4.077280e-06
## 7          healthy.000    acute_rejection.090 3.199889e-02 black 4.799833e-02
## 8          healthy.000    acute_rejection.365 6.882559e-01 black 6.882559e-01
## 9          healthy.000 acute_rejection.biopsy 3.350608e-01 black 3.518138e-01
## 10    no_rejection.000    acute_rejection.000 1.489642e-01 black 1.747311e-01
## 11    no_rejection.007    acute_rejection.007 2.290071e-02 black 4.007624e-02
## 12    no_rejection.090    acute_rejection.090 8.967787e-04 black 2.354044e-03
## 13    no_rejection.365    acute_rejection.365 1.013818e-03 black 2.365575e-03
## 16    no_rejection.090 acute_rejection.biopsy 5.687729e-03 black 1.085839e-02
## 17    no_rejection.365 acute_rejection.biopsy 2.736817e-02 black 4.421013e-02
## 18    no_rejection.000       no_rejection.007 1.152029e-05 blue3 3.456087e-05
## 21    no_rejection.007       no_rejection.090 3.593050e-02 blue3 5.030269e-02
## 23    no_rejection.090       no_rejection.365 1.755844e-01 blue3 1.940670e-01
## 24 acute_rejection.000    acute_rejection.007 1.390606e-03  red4 2.920273e-03
## 27 acute_rejection.007    acute_rejection.090 5.156805e-02  red4 6.768306e-02
## 29 acute_rejection.090    acute_rejection.365 1.497695e-01  red4 1.747311e-01

## only significant differences
pair_table %>% 
  filter(p<0.05) %>% 
  select(contrast, p, fdr) %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()
##                                           contrast            p          fdr
## 1                     healthy vs. no rejection D_0 1.507542e-15 3.165839e-14
## 2                     healthy vs. no rejection D_7 1.250083e-13 1.312587e-12
## 3                    healthy vs. no rejection D_90 5.134600e-09 2.695665e-08
## 4                   healthy vs. no rejection D_365 7.018490e-08 2.947766e-07
## 5                  healthy vs. acute rejection D_0 1.419560e-11 9.936919e-11
## 6                  healthy vs. acute rejection D_7 1.164937e-06 4.077280e-06
## 7                 healthy vs. acute rejection D_90 3.199889e-02 4.799833e-02
## 8         no rejection D_7 vs. acute rejection D_7 2.290071e-02 4.007624e-02
## 9       no rejection D_90 vs. acute rejection D_90 8.967787e-04 2.354044e-03
## 10    no rejection D_365 vs. acute rejection D_365 1.013818e-03 2.365575e-03
## 11  no rejection D_90 vs. acute rejection D_biopsy 5.687729e-03 1.085839e-02
## 12 no rejection D_365 vs. acute rejection D_biopsy 2.736817e-02 4.421013e-02
## 13           no rejection D_0 vs. no rejection D_7 1.152029e-05 3.456087e-05
## 14          no rejection D_7 vs. no rejection D_90 3.593050e-02 5.030269e-02
## 15     acute rejection D_0 vs. acute rejection D_7 1.390606e-03 2.920273e-03
##    star
## 1   ***
## 2   ***
## 3   ***
## 4   ***
## 5   ***
## 6   ***
## 7     *
## 8     *
## 9   ***
## 10   **
## 11   **
## 12    *
## 13  ***
## 14    *
## 15   **

pair_rest <- pair_table %>% 
  filter(p<0.05,
         term1 != 'healthy.000') %>% 
  mutate(star = if_else(p<0.05, '*', '')) %>% 
  mutate(star = if_else(p<0.01, '**', star)) %>% 
  mutate(star = if_else(p<0.001, '***', star)) %>%        
  data.frame()

Healthy does differ from all but acute rejection during biopsy and D365

3.7.3 Visualization

Open code

## plot
plotos <- data_long %>% 
  filter(!is.na(il36_b_value)) %>% 
  mutate(group = factor(group, levels = 
                          c("no_rejection", 
                            "acute_rejection",
                            "healthy")),
         time_point = factor(time_point)) %>% 
  mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>% 
  mutate(group_sep = interaction(group, time_point)) %>% 
  
## plotting
  ggplot(aes(x=time_point, y=il36_b_value)) +
  geom_boxplot(aes(fill=group, group = group_sep), outlier.shape = NA, 
               size=0.5,position = position_dodge(width = 0.84), color='black') +
  labs(x = "Time post TX", y = expression("IL-36" * beta * " (pg/ml)")) +
  scale_x_discrete(labels = c("Before TX", "7 days", "3 months", "365 days",
                              "Biopsy", "Healthy")) +
  theme(axis.text = element_text(size = 11),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 12)) +
  scale_fill_manual(values = c("no_rejection" = "skyblue2",
                               "acute_rejection" = "coral",
                               "healthy" = "darkolivegreen3"), 
                    labels = c("no rejection", "acute rejection", "healthy")) +
    coord_cartesian(ylim = c(0, 15.5)) +
  annotate("text", x = 6, y = 5, label = "#", size = 7) 

pair_rest$y <-  c(9.9,  7, 6.4,  9,  
                  8.1, 15.5, 12.5, 11.5)

plotos <- plotos + geom_segment(data = pair_rest, 
                                aes(x = xcoord_1, xend = xcoord_2, y = y, yend = y),
                                color = pair_rest$cole,
                                linetype = "solid", size = 0.4) 
plot_il36_b <- plotos + annotate("text", x = pair_rest$anot_coord, y = pair_rest$y+0.15, 
                            label = pair_rest$star, size = 7, color=pair_rest$cole) 
plot_il36_b

4 Correlation of ILs and renal function markers

4.1 Simple correlation matrix based on patient data

Open code
tr_cor <- patient_table %>%
   filter(group != 'healthy') %>% 
  select(egfr_mean:il36_b_mean) 


cor.mtest(tr_cor, conf.level = 0.95)
## $p
##                    egfr_mean creatinine_mean   il1_a_mean   il1_b_mean
## egfr_mean       0.000000e+00    2.103729e-16 6.139177e-01 2.941143e-01
## creatinine_mean 2.103729e-16    0.000000e+00 9.054514e-01 9.825541e-01
## il1_a_mean      6.139177e-01    9.054514e-01 0.000000e+00 5.817931e-18
## il1_b_mean      2.941143e-01    9.825541e-01 5.817931e-18 0.000000e+00
## il1_ra_mean     7.275104e-01    8.430982e-01 4.288839e-15 2.316736e-11
## il18_mean       8.565717e-01    7.798943e-01 1.680311e-15 7.659848e-14
## il18_bp_mean    1.795693e-01    2.157285e-01 3.841527e-01 7.085025e-01
## il18_free_mean  6.704398e-01    6.534084e-01 3.949705e-15 1.055580e-11
## il36_b_mean     2.452979e-01    3.918567e-01 9.302473e-05 8.379317e-06
##                  il1_ra_mean    il18_mean il18_bp_mean il18_free_mean
## egfr_mean       7.275104e-01 8.565717e-01    0.1795693   6.704398e-01
## creatinine_mean 8.430982e-01 7.798943e-01    0.2157285   6.534084e-01
## il1_a_mean      4.288839e-15 1.680311e-15    0.3841527   3.949705e-15
## il1_b_mean      2.316736e-11 7.659848e-14    0.7085025   1.055580e-11
## il1_ra_mean     0.000000e+00 2.176954e-20    0.5645638   2.259573e-13
## il18_mean       2.176954e-20 0.000000e+00    0.9768424   6.833531e-85
## il18_bp_mean    5.645638e-01 9.768424e-01    0.0000000   4.524752e-01
## il18_free_mean  2.259573e-13 6.833531e-85    0.4524752   0.000000e+00
## il36_b_mean     2.540206e-02 1.512619e-03    0.3023744   1.849524e-03
##                  il36_b_mean
## egfr_mean       2.452979e-01
## creatinine_mean 3.918567e-01
## il1_a_mean      9.302473e-05
## il1_b_mean      8.379317e-06
## il1_ra_mean     2.540206e-02
## il18_mean       1.512619e-03
## il18_bp_mean    3.023744e-01
## il18_free_mean  1.849524e-03
## il36_b_mean     0.000000e+00
## 
## $lowCI
##                   egfr_mean creatinine_mean il1_a_mean  il1_b_mean il1_ra_mean
## egfr_mean        1.00000000     -0.71827559 -0.1246886 -0.07833667 -0.19605109
## creatinine_mean -0.71827559      1.00000000 -0.1770066 -0.16527820 -0.18358576
## il1_a_mean      -0.12468860     -0.17700661  1.0000000  0.54243007  0.48621354
## il1_b_mean      -0.07833667     -0.16527820  0.5424301  1.00000000  0.39813645
## il1_ra_mean     -0.19605109     -0.18358576  0.4862135  0.39813645  1.00000000
## il18_mean       -0.18215813     -0.14367730  0.4947200  0.45880699  0.58405262
## il18_bp_mean    -0.27670352     -0.06215646 -0.2387838 -0.19815124 -0.21478835
## il18_free_mean  -0.20241378     -0.12937782  0.4869689  0.40705248  0.44800513
## il36_b_mean     -0.06868499     -0.23764656  0.1685782  0.21530828  0.02392036
##                  il18_mean il18_bp_mean il18_free_mean il36_b_mean
## egfr_mean       -0.1821581  -0.27670352     -0.2024138 -0.06868499
## creatinine_mean -0.1436773  -0.06215646     -0.1293778 -0.23764656
## il1_a_mean       0.4947200  -0.23878380      0.4869689  0.16857815
## il1_b_mean       0.4588070  -0.19815124      0.4070525  0.21530828
## il1_ra_mean      0.5840526  -0.21478835      0.4480051  0.02392036
## il18_mean        1.0000000  -0.16467977      0.9574656  0.10512859
## il18_bp_mean    -0.1646798   1.00000000     -0.2291107 -0.07985789
## il18_free_mean   0.9574656  -0.22911074      1.0000000  0.10003378
## il36_b_mean      0.1051286  -0.07985789      0.1000338  1.00000000
## 
## $uppCI
##                   egfr_mean creatinine_mean il1_a_mean il1_b_mean il1_ra_mean
## egfr_mean        1.00000000     -0.51293155 0.20891140  0.2532441   0.1378676
## creatinine_mean -0.51293155      1.00000000 0.15716924  0.1689303   0.1505300
## il1_a_mean       0.20891140      0.15716924 1.00000000  0.7374956   0.7005925
## il1_b_mean       0.25324409      0.16893033 0.73749557  1.0000000   0.6403732
## il1_ra_mean      0.13786762      0.15052999 0.70059252  0.6403732   1.0000000
## il18_mean        0.15197323      0.19034518 0.70625117  0.6821760   0.7640880
## il18_bp_mean     0.05321464      0.26839831 0.09361479  0.1357235   0.1186266
## il18_free_mean   0.13136199      0.20434877 0.70109612  0.6466077   0.6748387
## il36_b_mean      0.26230338      0.09480974 0.46805481  0.5050816   0.3463574
##                 il18_mean il18_bp_mean il18_free_mean il36_b_mean
## egfr_mean       0.1519732   0.05321464      0.1313620  0.26230338
## creatinine_mean 0.1903452   0.26839831      0.2043488  0.09480974
## il1_a_mean      0.7062512   0.09361479      0.7010961  0.46805481
## il1_b_mean      0.6821760   0.13572347      0.6466077  0.50508160
## il1_ra_mean     0.7640880   0.11862657      0.6748387  0.34635737
## il18_mean       1.0000000   0.16952788      0.9781094  0.41603888
## il18_bp_mean    0.1695279   1.00000000      0.1037483  0.25181091
## il18_free_mean  0.9781094   0.10374828      1.0000000  0.41177200
## il36_b_mean     0.4160389   0.25181091      0.4117720  1.00000000

tr_cor <- cor(tr_cor, method='spearman', 
              use = 'pairwise.complete.obs') 
corrplot(tr_cor, addCoef.col = 'black',tl.col = "white")
text(1:9, 9.9, expression("eGFR", "Creatinine", "IL-1" * alpha, 
                           "IL-1" * beta, "IL-1 RA", "IL-18", "IL-18 BP",
               "IL-18 (free)",  "IL-36" * beta))

text(-0.1, 9:1,  expression("eGFR", "Creatinine", "IL-1" * alpha, 
                           "IL-1" * beta, "IL-1 RA", "IL-18", "IL-18 BP",
               "IL-18 (free)",  "IL-1" * beta))

There is not any significant correlation between interleukins and markers of renal function. At least modest (but still not statistically significant) correlation is for interleukin-18BP.

4.2 Adjsuting IL-18BP model for renal function

To further validate that renal functions does not play an important role for above-shown findings, we will re-fit mixed-effects model adjusting also for the effect of renal function. We will do it only for IL-18BP where modest correlation between average values of renal markers and IL-18BP levels was found.

4.2.1 Adjustment for EGFR

Open code
model_data <- data_long %>% 
  filter(!is.na(il18_bp_value))

## Complex models (incuding healthy and biopsy data)
model_il18_bp <-  lmer(log10(il18_bp_value) ~ 
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_bp_adjEGFR <-  lmer(log10(il18_bp_value) ~ log2(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_bp_adjEGFR_REML <-  lmer(log10(il18_bp_value) ~ log2(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_bp_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_bp_value)
##                 F Df Df.res   Pr(>F)   
## log2(egfr) 7.9151  1 547.11 0.005079 **
## group_sep  0.8430  9 498.72 0.576683   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il18_bp_adjEGFR_REML)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_bp_value) ~ log2(egfr) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: -484.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.0627  -0.3090   0.0389   0.3262   4.7340 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.00581  0.07622 
##  Residual             0.01870  0.13675 
## Number of obs: 577, groups:  id, 169
## 
## Fixed effects:
##                                  Estimate Std. Error t value
## (Intercept)                      3.541856   0.035691  99.236
## log2(egfr)                      -0.023511   0.008321  -2.825
## group_sephealthy.000            -0.082453   0.047640  -1.731
## group_sepno_rejection.000        0.003293   0.029488   0.112
## group_sepacute_rejection.007    -0.022998   0.032948  -0.698
## group_sepno_rejection.007       -0.008345   0.032425  -0.257
## group_sepacute_rejection.090    -0.061574   0.037176  -1.656
## group_sepno_rejection.090       -0.043781   0.037677  -1.162
## group_sepacute_rejection.365    -0.030476   0.038551  -0.791
## group_sepno_rejection.365       -0.054496   0.038220  -1.426
## group_sepacute_rejection.biopsy -0.015365   0.036709  -0.419
## 
## Correlation of Fixed Effects:
##              (Intr) lg2(g) g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log2(egfr)    0.724                                                
## grp_sph.000  -0.805 -0.620                                         
## grp_sp_.000  -0.568  0.011  0.425                                  
## grp_spc_.007 -0.703 -0.414  0.559  0.484                           
## grp_spn_.007 -0.827 -0.419  0.652  0.694   0.617                   
## grp_spc_.090 -0.769 -0.568  0.620  0.427   0.623        0.632      
## grp_spn_.090 -0.895 -0.613  0.718  0.595   0.636        0.804      
## grp_spc_.365 -0.755 -0.566  0.610  0.412   0.609        0.618      
## grp_spn_.365 -0.894 -0.621  0.718  0.586   0.634        0.799      
## grp_spct_r.  -0.622 -0.364  0.494  0.430   0.539        0.547      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log2(egfr)                                                      
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.687                                             
## grp_spc_.365  0.653        0.675                                
## grp_spn_.365  0.687        0.845        0.674                   
## grp_spct_r.   0.551        0.563        0.535        0.561
AIC(model_il18_bp, model_il18_bp_adjEGFR)
##                       df       AIC
## model_il18_bp         12 -519.0147
## model_il18_bp_adjEGFR 13 -525.0759
BIC(model_il18_bp, model_il18_bp_adjEGFR)
##                       df       BIC
## model_il18_bp         12 -466.7206
## model_il18_bp_adjEGFR 13 -468.4240


## Patient-only, planned-measures-only models
model_data <- data_long %>% 
  filter(!is.na(il18_bp_value),
         group != 'healthy',
         time_point != 'biopsy')

model_il18_bp <-  lmer(log10(il18_bp_value) ~ 
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_bp_adjEGFR <-  lmer(log10(il18_bp_value) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_bp_adjEGFR_REML <-  lmer(log10(il18_bp_value) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_bp_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_bp_value)
##                       F Df Df.res   Pr(>F)   
## log2(egfr)       8.3333  1 493.34 0.004063 **
## group            0.0399  1 136.84 0.841907   
## time_point       1.1440  3 427.71 0.330982   
## group:time_point 0.4743  3 383.95 0.700327   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il18_bp_adjEGFR_REML)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_bp_value) ~ log2(egfr) + group * time_point + (1 |      id)
##    Data: model_data
## 
## REML criterion at convergence: -395.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.6489 -0.3011  0.0298  0.3526  4.5301 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.006333 0.07958 
##  Residual             0.020298 0.14247 
## Number of obs: 518, groups:  id, 138
## 
## Fixed effects:
##                                     Estimate Std. Error t value
## (Intercept)                         3.537355   0.032864 107.637
## log2(egfr)                         -0.025990   0.008963  -2.900
## groupacute_rejection               -0.003096   0.030741  -0.101
## time_point007                      -0.007498   0.025666  -0.292
## time_point090                      -0.040094   0.033041  -1.213
## time_point365                      -0.050629   0.033752  -1.500
## groupacute_rejection:time_point007 -0.011536   0.037562  -0.307
## groupacute_rejection:time_point090 -0.015374   0.038446  -0.400
## groupacute_rejection:time_point365  0.027516   0.039657   0.694
## 
## Correlation of Fixed Effects:
##             (Intr) lg2(g) grpct_ tm_007 tm_090 tm_365 g_:_00 g_:_09
## log2(egfr)   0.857                                                 
## grpct_rjctn -0.293 -0.012                                          
## time_pnt007 -0.759 -0.583  0.284                                   
## time_pnt090 -0.856 -0.764  0.224  0.703                            
## time_pnt365 -0.855 -0.768  0.220  0.700  0.783                     
## grpct_:_007  0.184  0.008 -0.635 -0.455 -0.182 -0.178              
## grpct_:_090  0.230  0.066 -0.622 -0.260 -0.409 -0.220  0.509       
## grpct_:_365  0.218  0.058 -0.604 -0.249 -0.212 -0.395  0.494  0.486
AIC(model_il18_bp, model_il18_bp_adjEGFR)
##                       df       AIC
## model_il18_bp         10 -422.7043
## model_il18_bp_adjEGFR 11 -429.1605
BIC(model_il18_bp, model_il18_bp_adjEGFR)
##                       df       BIC
## model_il18_bp         10 -380.2045
## model_il18_bp_adjEGFR 11 -382.4108

Inclusion of EGFR improves the fit noticeably

4.2.2 Adjustment for creatinine

Open code


## Complex models (including healthy and biopsy data)

model_data <- data_long %>% 
  filter(!is.na(il18_bp_value))

model_il18_bp <-  lmer(log10(il18_bp_value) ~ 
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_bp_adjCRE <-  lmer(log2(il18_bp_value) ~ log2(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_bp_adjCRE_REML <-  lmer(log10(il18_bp_value) ~ log2(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_bp_adjCRE_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_bp_value)
##                       F Df Df.res   Pr(>F)   
## log2(creatinine) 8.4351  1 534.24 0.003833 **
## group_sep        0.8406  9 499.16 0.578871   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il18_bp_adjCRE_REML)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il18_bp_value) ~ log2(creatinine) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: -485.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.0607  -0.3099   0.0346   0.3382   4.7544 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.005777 0.0760  
##  Residual             0.018700 0.1367  
## Number of obs: 577, groups:  id, 169
## 
## Fixed effects:
##                                  Estimate Std. Error t value
## (Intercept)                      3.345720   0.095493  35.036
## log2(creatinine)                 0.028785   0.009868   2.917
## group_sephealthy.000            -0.080801   0.047398  -1.705
## group_sepno_rejection.000        0.004245   0.029466   0.144
## group_sepacute_rejection.007    -0.022510   0.032840  -0.685
## group_sepno_rejection.007       -0.006920   0.032434  -0.213
## group_sepacute_rejection.090    -0.060569   0.036995  -1.637
## group_sepno_rejection.090       -0.041914   0.037610  -1.114
## group_sepacute_rejection.365    -0.029155   0.038422  -0.759
## group_sepno_rejection.365       -0.052324   0.038205  -1.370
## group_sepacute_rejection.biopsy -0.014570   0.036653  -0.398
## 
## Correlation of Fixed Effects:
##              (Intr) lg2(c) g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log2(crtnn)  -0.966                                                
## grp_sph.000  -0.729  0.616                                         
## grp_sp_.000  -0.216  0.001  0.434                                  
## grp_spc_.007 -0.545  0.408  0.556  0.490                           
## grp_spn_.007 -0.602  0.421  0.653  0.698   0.617                   
## grp_spc_.090 -0.678  0.562  0.617  0.436   0.620        0.632      
## grp_spn_.090 -0.760  0.612  0.716  0.602   0.633        0.804      
## grp_spc_.365 -0.672  0.562  0.607  0.421   0.606        0.618      
## grp_spn_.365 -0.766  0.621  0.717  0.593   0.631        0.799      
## grp_spct_r.  -0.483  0.360  0.492  0.435   0.537        0.547      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log2(crtnn)                                                     
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.685                                             
## grp_spc_.365  0.650        0.673                                
## grp_spn_.365  0.685        0.845        0.673                   
## grp_spct_r.   0.549        0.561        0.533        0.559
AIC(model_il18_bp, model_il18_bp_adjCRE)
##                      df       AIC
## model_il18_bp        12 -519.0147
## model_il18_bp_adjCRE 13  859.8325
BIC(model_il18_bp, model_il18_bp_adjCRE)
##                      df       BIC
## model_il18_bp        12 -466.7206
## model_il18_bp_adjCRE 13  916.4844


## Patient-only, planned-measures-only models
model_data <- data_long %>% 
  filter(!is.na(il18_bp_value),
         group != 'healthy',
         time_point != 'biopsy')

model_il18_bp <-  lmer(log10(il18_bp_value) ~ 
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_bp_adjEGFR <-  lmer(log10(il18_bp_value) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_bp_adjEGFR_REML <-  lmer(log10(il18_bp_value) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_bp_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_bp_value)
##                       F Df Df.res   Pr(>F)   
## log2(egfr)       8.3333  1 493.34 0.004063 **
## group            0.0399  1 136.84 0.841907   
## time_point       1.1440  3 427.71 0.330982   
## group:time_point 0.4743  3 383.95 0.700327   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_bp, model_il18_bp_adjEGFR)
##                       df       AIC
## model_il18_bp         10 -422.7043
## model_il18_bp_adjEGFR 11 -429.1605
BIC(model_il18_bp, model_il18_bp_adjEGFR)
##                       df       BIC
## model_il18_bp         10 -380.2045
## model_il18_bp_adjEGFR 11 -382.4108

The analysis shows that when renal function are accounted for (included as a covariate to the mixed-effects model), there is no remaining variability explainable by other factors such as group [acute rejection/no rejection/healthy], time_point and group:timepoint interaction. It suggests that IL-18BP might be driven mainly with renal function and its change over time whereas the effects of group and timepoint etc may be indirect.

4.3 Adjusted free IL-18 model

Given that IL-18BP is affected with renal function, and the IL18-BP affects level of other forms of IL-18 (including free IL-18), we will fit model explaining the IL-18 levels but simultaneously adjusting for renal functions

4.3.1 Adjustment for EGFR

Open code
model_data <- data_long %>% 
  filter(!is.na(il18_free))

## Complex models (incuding healthy and biopsy data)
model_il18_free <-  lmer(log10(il18_free) ~ 
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_free_adjEGFR <-  lmer(log10(il18_free) ~ log2(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_free_adjEGFR_REML <-  lmer(log10(il18_free) ~ log2(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_free_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_free)
##                  F Df Df.res Pr(>F)    
## log2(egfr)  0.6484  1  454.0 0.4211    
## group_sep  18.2407  9  408.3 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_free, model_il18_free_adjEGFR)
##                         df      AIC
## model_il18_free         12 227.9601
## model_il18_free_adjEGFR 13 229.2903
BIC(model_il18_free, model_il18_free_adjEGFR)
##                         df      BIC
## model_il18_free         12 278.0706
## model_il18_free_adjEGFR 13 283.5766

## contrasts
em_il18_free <- emmeans(model_il18_free_adjEGFR_REML, specs = pairwise ~ group_sep, 
              adjust = 'none', type = 'response')

em_il18_free_cont <- data.frame(em_il18_free$contrasts)[c(
  1:3, 10:17, 19, 25, 26, 32, 36, 37, 39, 42, 43, 45
),]

as_tibble(em_il18_free_cont) %>% 
  filter(p.value<0.05) %>% data.frame()
##                                     contrast     ratio         SE       df null
## 1          acute_rejection.000 / healthy.000 2.4073392 0.53767731 434.9662    1
## 2     acute_rejection.000 / no_rejection.000 0.6793931 0.09351587 425.4130    1
## 3  acute_rejection.000 / acute_rejection.007 1.7959185 0.27628214 368.5469    1
## 4             healthy.000 / no_rejection.000 0.2822174 0.05813244 432.7976    1
## 5             healthy.000 / no_rejection.007 0.6344966 0.10611335 418.2433    1
## 6          healthy.000 / acute_rejection.090 0.3025239 0.05310212 416.0604    1
## 7             healthy.000 / no_rejection.090 0.3225245 0.05607866 443.6156    1
## 8          healthy.000 / acute_rejection.365 0.3762615 0.06780092 424.0318    1
## 9             healthy.000 / no_rejection.365 0.4349602 0.07537710 443.7478    1
## 10      healthy.000 / acute_rejection.biopsy 0.3113010 0.06239168 447.3936    1
## 11       no_rejection.000 / no_rejection.007 2.2482543 0.25255124 406.6198    1
## 12 acute_rejection.007 / acute_rejection.090 0.4055182 0.05689067 337.8121    1
## 13       no_rejection.007 / no_rejection.090 0.5083156 0.06651835 391.4777    1
##      t.ratio      p.value
## 1   3.933401 9.745268e-05
## 2  -2.808326 5.209278e-03
## 3   3.806037 1.653100e-04
## 4  -6.141612 1.851593e-09
## 5  -2.720179 6.796665e-03
## 6  -6.811330 3.398621e-11
## 7  -6.508020 2.061651e-10
## 8  -5.424479 9.787394e-08
## 9  -4.803908 2.132018e-06
## 10 -5.822679 1.106736e-08
## 11  7.212129 2.710851e-12
## 12 -6.433682 4.262068e-10
## 13 -5.170801 3.726822e-07

## Patient-only, planned-measures-only models
model_data <- data_long %>% 
  filter(!is.na(il18_free),
         group != 'healthy',
         time_point != 'biopsy')

model_il18_free <-  lmer(log10(il18_free) ~ 
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_free_adjEGFR <-  lmer(log10(il18_free) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_free_adjEGFR_REML <-  lmer(log10(il18_free) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_free_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_free)
##                        F Df Df.res  Pr(>F)    
## log2(egfr)        1.1034  1 401.46 0.29415    
## group             1.7354  1 127.35 0.19008    
## time_point       38.0225  3 333.84 < 2e-16 ***
## group:time_point  2.9089  3 300.10 0.03485 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_free, model_il18_free_adjEGFR)
##                         df      AIC
## model_il18_free         10 228.2280
## model_il18_free_adjEGFR 11 229.0902
BIC(model_il18_free, model_il18_free_adjEGFR)
##                         df      BIC
## model_il18_free         10 268.6781
## model_il18_free_adjEGFR 11 273.5853

Inclusion of EGFR improves the fit noticeably

4.3.2 Adjustment for creatinine

Open code


## Complex models (including healthy and biopsy data)

model_data <- data_long %>% 
  filter(!is.na(il18_free))

model_il18_free <-  lmer(log10(il18_free) ~ 
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_free_adjCRE <-  lmer(log10(il18_free) ~ log2(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_free_adjCRE_REML <-  lmer(log10(il18_free) ~ log2(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_free_adjCRE_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_free)
##                        F Df Df.res Pr(>F)    
## log2(creatinine)  0.0606  1 445.85 0.8057    
## group_sep        18.0057  9 408.47 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_free, model_il18_free_adjCRE)
##                        df      AIC
## model_il18_free        12 227.9601
## model_il18_free_adjCRE 13 229.8973
BIC(model_il18_free, model_il18_free_adjCRE)
##                        df      BIC
## model_il18_free        12 278.0706
## model_il18_free_adjCRE 13 284.1836

## contrasts
em_il18_free <- emmeans(model_il18_free_adjCRE_REML, specs = pairwise ~ group_sep, 
              adjust = 'none', type = 'response')

em_il18_free_cont <- data.frame(em_il18_free$contrasts)[c(
  1:3, 10:17, 19, 25, 26, 32, 36, 37, 39, 42, 43, 45
),]

as_tibble(em_il18_free_cont) %>% 
  filter(p.value<0.05) %>% data.frame()
##                                     contrast     ratio         SE       df null
## 1          acute_rejection.000 / healthy.000 2.2241926 0.49446636 431.4742    1
## 2     acute_rejection.000 / no_rejection.000 0.6801853 0.09368936 425.2481    1
## 3  acute_rejection.000 / acute_rejection.007 1.7300734 0.26527877 368.8704    1
## 4             healthy.000 / no_rejection.000 0.3058122 0.06224851 427.1593    1
## 5             healthy.000 / no_rejection.007 0.6623566 0.11018206 415.2399    1
## 6          healthy.000 / acute_rejection.090 0.3093849 0.05432380 414.9301    1
## 7             healthy.000 / no_rejection.090 0.3284553 0.05704034 442.9723    1
## 8          healthy.000 / acute_rejection.365 0.3842716 0.06922734 423.0664    1
## 9             healthy.000 / no_rejection.365 0.4421845 0.07652278 443.2303    1
## 10      healthy.000 / acute_rejection.biopsy 0.3250711 0.06504205 445.2994    1
## 11       no_rejection.000 / no_rejection.007 2.1658930 0.24151755 408.1660    1
## 12 acute_rejection.007 / acute_rejection.090 0.3977471 0.05580360 337.5963    1
## 13       no_rejection.007 / no_rejection.090 0.4958889 0.06485914 391.3846    1
##      t.ratio      p.value
## 1   3.595808 3.607667e-04
## 2  -2.797933 5.376774e-03
## 3   3.574970 3.968385e-04
## 4  -5.820564 1.152945e-08
## 5  -2.476434 1.366773e-02
## 6  -6.681433 7.635869e-11
## 7  -6.411028 3.708165e-10
## 8  -5.308879 1.784928e-07
## 9  -4.715392 3.236763e-06
## 10 -5.616153 3.443165e-08
## 11  6.930648 1.642603e-11
## 12 -6.571234 1.892147e-10
## 13 -5.362670 1.405310e-07


## Patient-only, planned-measures-only models
model_data <- data_long %>% 
  filter(!is.na(il18_free),
         group != 'healthy',
         time_point != 'biopsy')

model_il18_free <-  lmer(log10(il18_free) ~ 
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_free_adjEGFR <-  lmer(log10(il18_free) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_free_adjEGFR_REML <-  lmer(log10(il18_free) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_free_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_free)
##                        F Df Df.res  Pr(>F)    
## log2(egfr)        1.1034  1 401.46 0.29415    
## group             1.7354  1 127.35 0.19008    
## time_point       38.0225  3 333.84 < 2e-16 ***
## group:time_point  2.9089  3 300.10 0.03485 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_free, model_il18_free_adjEGFR)
##                         df      AIC
## model_il18_free         10 228.2280
## model_il18_free_adjEGFR 11 229.0902
BIC(model_il18_free, model_il18_free_adjEGFR)
##                         df      BIC
## model_il18_free         10 268.6781
## model_il18_free_adjEGFR 11 273.5853

Adjustment for renal function does not modify result of the model explaining free IL-18. Markers of renal function do not clearly affect the pattern of difference in free IL-18 levels across groups or time points

What about total IL-18?

4.4 Adjusted total IL-18 model

4.4.1 Adjustment for EGFR

Open code
model_data <- data_long %>% 
  filter(!is.na(il18_value))

## Complex models (incuding healthy and biopsy data)
model_il18_value <-  lmer(log10(il18_value) ~ 
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_value_adjEGFR <-  lmer(log10(il18_value) ~ log2(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_value_adjEGFR_REML <-  lmer(log10(il18_value) ~ log2(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_value_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_value)
##                  F Df Df.res Pr(>F)    
## log2(egfr)  0.0518  1 479.60 0.8201    
## group_sep  20.2839  9 420.04 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_value, model_il18_value_adjEGFR)
##                          df      AIC
## model_il18_value         12 271.2115
## model_il18_value_adjEGFR 13 273.1579
BIC(model_il18_value, model_il18_value_adjEGFR)
##                          df      BIC
## model_il18_value         12 321.7868
## model_il18_value_adjEGFR 13 327.9478

## contrasts
em_il18_value <- emmeans(model_il18_value_adjEGFR_REML, specs = pairwise ~ group_sep, 
              adjust = 'none', type = 'response')

em_il18_value_cont <- data.frame(em_il18_value$contrasts)[c(
  1:3, 10:17, 19, 25, 26, 32, 36, 37, 39, 42, 43, 45
),]

as_tibble(em_il18_value_cont) %>% 
  filter(p.value<0.05) %>% data.frame()
##                                     contrast     ratio         SE       df null
## 1          acute_rejection.000 / healthy.000 2.8016661 0.60639470 454.4073    1
## 2     acute_rejection.000 / no_rejection.000 0.6587990 0.09349772 431.3921    1
## 3  acute_rejection.000 / acute_rejection.007 1.8091393 0.28340410 375.4448    1
## 4             healthy.000 / no_rejection.000 0.2351455 0.04639809 455.5534    1
## 5          healthy.000 / acute_rejection.007 0.6457369 0.11290118 430.7954    1
## 6             healthy.000 / no_rejection.007 0.5400838 0.08320735 434.4828    1
## 7          healthy.000 / acute_rejection.090 0.2686464 0.04407996 427.5350    1
## 8             healthy.000 / no_rejection.090 0.2750560 0.04449238 467.0760    1
## 9          healthy.000 / acute_rejection.365 0.3259036 0.05519896 439.1190    1
## 10            healthy.000 / no_rejection.365 0.3689904 0.05944705 467.1907    1
## 11      healthy.000 / acute_rejection.biopsy 0.2703141 0.05182941 468.8662    1
## 12       no_rejection.000 / no_rejection.007 2.2968076 0.26377620 414.3946    1
## 13 acute_rejection.007 / acute_rejection.090 0.4160308 0.05972821 342.3573    1
## 14       no_rejection.007 / no_rejection.090 0.5092838 0.06839126 396.4367    1
##      t.ratio      p.value
## 1   4.759798 2.608696e-06
## 2  -2.940618 3.451797e-03
## 3   3.784527 1.791956e-04
## 4  -7.336186 1.014798e-12
## 5  -2.501493 1.273714e-02
## 6  -3.998545 7.488620e-05
## 7  -8.010396 1.096085e-14
## 8  -7.979724 1.144495e-14
## 9  -6.619473 1.053554e-10
## 10 -6.188327 1.330858e-09
## 11 -6.822709 2.768522e-11
## 12  7.240387 2.192254e-12
## 13 -6.108626 2.730562e-09
## 14 -5.024607 7.647530e-07

## Patient-only, planned-measures-only models
model_data <- data_long %>% 
  filter(!is.na(il18_value),
         group != 'healthy',
         time_point != 'biopsy')

model_il18_value <-  lmer(log10(il18_value) ~ 
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_value_adjEGFR <-  lmer(log10(il18_value) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_value_adjEGFR_REML <-  lmer(log10(il18_value) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_value_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_value)
##                        F Df Df.res  Pr(>F)    
## log2(egfr)        0.3786  1 404.44 0.53871    
## group             2.4964  1 127.58 0.11658    
## time_point       38.6155  3 335.53 < 2e-16 ***
## group:time_point  3.0134  3 301.04 0.03035 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_value, model_il18_value_adjEGFR)
##                          df      AIC
## model_il18_value         10 229.1316
## model_il18_value_adjEGFR 11 230.7408
BIC(model_il18_value, model_il18_value_adjEGFR)
##                          df      BIC
## model_il18_value         10 269.6289
## model_il18_value_adjEGFR 11 275.2879

4.4.2 Adjustment for creatinine

Open code


## Complex models (including healthy and biopsy data)

model_data <- data_long %>% 
  filter(!is.na(il18_value))

model_il18_value <-  lmer(log10(il18_value) ~ 
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_value_adjCRE <-  lmer(log10(il18_value) ~ log2(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = FALSE)

model_il18_value_adjCRE_REML <-  lmer(log10(il18_value) ~ log2(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_value_adjCRE_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_value)
##                        F Df Df.res Pr(>F)    
## log2(creatinine)  0.0663  1 471.87 0.7969    
## group_sep        20.0253  9 420.33 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_value, model_il18_value_adjCRE)
##                         df      AIC
## model_il18_value        12 271.2115
## model_il18_value_adjCRE 13 273.1441
BIC(model_il18_value, model_il18_value_adjCRE)
##                         df      BIC
## model_il18_value        12 321.7868
## model_il18_value_adjCRE 13 327.9341

## contrasts
em_il18_value <- emmeans(model_il18_value_adjCRE_REML, specs = pairwise ~ group_sep, 
              adjust = 'none', type = 'response')

em_il18_value_cont <- data.frame(em_il18_value$contrasts)[c(
  1:3, 10:17, 19, 25, 26, 32, 36, 37, 39, 42, 43, 45
),]

as_tibble(em_il18_value_cont) %>% 
  filter(p.value<0.05) %>% data.frame()
##                                     contrast     ratio         SE       df null
## 1          acute_rejection.000 / healthy.000 2.6120988 0.56418523 450.2078    1
## 2     acute_rejection.000 / no_rejection.000 0.6587692 0.09347288 431.7167    1
## 3  acute_rejection.000 / acute_rejection.007 1.7506807 0.27327263 376.1026    1
## 4             healthy.000 / no_rejection.000 0.2521992 0.04927002 448.6505    1
## 5          healthy.000 / acute_rejection.007 0.6702199 0.11738276 429.1309    1
## 6             healthy.000 / no_rejection.007 0.5606891 0.08608655 430.9211    1
## 7          healthy.000 / acute_rejection.090 0.2740685 0.04503086 426.8454    1
## 8             healthy.000 / no_rejection.090 0.2794271 0.04516298 466.4659    1
## 9          healthy.000 / acute_rejection.365 0.3320203 0.05626098 438.3287    1
## 10            healthy.000 / no_rejection.365 0.3741927 0.06021292 466.7281    1
## 11      healthy.000 / acute_rejection.biopsy 0.2808052 0.05384462 466.7466    1
## 12       no_rejection.000 / no_rejection.007 2.2231996 0.25325311 416.6459    1
## 13 acute_rejection.007 / acute_rejection.090 0.4089233 0.05870204 342.3214    1
## 14       no_rejection.007 / no_rejection.090 0.4983638 0.06686930 396.7331    1
##      t.ratio      p.value
## 1   4.445379 1.106241e-05
## 2  -2.941585 3.441110e-03
## 3   3.587587 3.778455e-04
## 4  -7.051215 6.743530e-12
## 5  -2.284731 2.281589e-02
## 6  -3.768398 1.871701e-04
## 7  -7.877888 2.790387e-14
## 8  -7.888616 2.190487e-14
## 9  -6.506677 2.102148e-10
## 10 -6.108748 2.118180e-09
## 11 -6.623671 9.669253e-11
## 12  7.013614 9.438174e-12
## 13 -6.229265 1.374605e-09
## 14 -5.190319 3.358860e-07


## Patient-only, planned-measures-only models
model_data <- data_long %>% 
  filter(!is.na(il18_value),
         group != 'healthy',
         time_point != 'biopsy')

model_il18_value <-  lmer(log10(il18_value) ~ 
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_value_adjEGFR <-  lmer(log10(il18_value) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = FALSE)

model_il18_value_adjEGFR_REML <-  lmer(log10(il18_value) ~ log2(egfr) +
                       group*time_point +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il18_value_adjEGFR_REML, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il18_value)
##                        F Df Df.res  Pr(>F)    
## log2(egfr)        0.3786  1 404.44 0.53871    
## group             2.4964  1 127.58 0.11658    
## time_point       38.6155  3 335.53 < 2e-16 ***
## group:time_point  3.0134  3 301.04 0.03035 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_il18_value, model_il18_value_adjEGFR)
##                          df      AIC
## model_il18_value         10 229.1316
## model_il18_value_adjEGFR 11 230.7408
BIC(model_il18_value, model_il18_value_adjEGFR)
##                          df      BIC
## model_il18_value         10 269.6289
## model_il18_value_adjEGFR 11 275.2879

Adjustment for renal function does not modify the results of the models explaining total IL-18. Markers of renal function do not clearly affect the above-shown patterns in total IL-18 levels and their difference across groups or time points

4.5 Renal function adjusment for other interleukins

4.5.1 IL-1a

Open code
model_data <- data_long %>% 
  filter(!is.na(il1_a_value))

## eGFR
model_il1a_value_adjEGFR <-  lmer(log10(il1_a_value) ~ log10(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il1a_value_adjEGFR, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_a_value)
##                   F Df Df.res Pr(>F)    
## log10(egfr)  0.3496  1 447.83 0.5547    
## group_sep   14.1460  9 424.92 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_a_value) ~ log10(egfr) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 810.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5716 -0.3349  0.0666  0.5774  3.5006 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.05197  0.228   
##  Residual             0.23910  0.489   
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      0.59923    0.12586   4.761
## log10(egfr)                      0.05827    0.09805   0.594
## group_sephealthy.000            -1.02366    0.15615  -6.556
## group_sepno_rejection.000        0.01424    0.10245   0.139
## group_sepacute_rejection.007    -0.11283    0.11920  -0.947
## group_sepno_rejection.007       -0.12177    0.11330  -1.075
## group_sepacute_rejection.090    -0.06818    0.13322  -0.512
## group_sepno_rejection.090        0.14640    0.14589   1.004
## group_sepacute_rejection.365    -0.16452    0.13886  -1.185
## group_sepno_rejection.365        0.08185    0.14728   0.556
## group_sepacute_rejection.biopsy -0.08932    0.13175  -0.678
## 
## Correlation of Fixed Effects:
##              (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log10(egfr)   0.730                                                
## grp_sph.000  -0.865 -0.669                                         
## grp_sp_.000  -0.569  0.007  0.458                                  
## grp_spc_.007 -0.716 -0.415  0.611  0.505                           
## grp_spn_.007 -0.831 -0.427  0.704  0.682   0.636                   
## grp_spc_.090 -0.783 -0.563  0.676  0.452   0.625        0.653      
## grp_spn_.090 -0.811 -0.559  0.699  0.528   0.589        0.720      
## grp_spc_.365 -0.765 -0.559  0.662  0.435   0.608        0.635      
## grp_spn_.365 -0.816 -0.571  0.704  0.523   0.590        0.720      
## grp_spct_r.  -0.637 -0.359  0.542  0.458   0.543        0.569      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log10(egfr)                                                     
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.635                                             
## grp_spc_.365  0.652        0.621                                
## grp_spn_.365  0.639        0.715        0.624                   
## grp_spct_r.   0.556        0.524        0.537        0.525

## creatinine
model_il1a_value_adjCRE <-  lmer(log10(il1_a_value) ~ log10(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il1a_value_adjCRE, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_a_value)
##                         F Df Df.res Pr(>F)    
## log10(creatinine)  0.2582  1 438.05 0.6116    
## group_sep         14.1083  9 424.95 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjCRE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_a_value) ~ log10(creatinine) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 810
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5825 -0.3292  0.0701  0.5676  3.5070 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.05217  0.2284  
##  Residual             0.23901  0.4889  
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      0.71098    0.33694   2.110
## log10(creatinine)               -0.05901    0.11556  -0.511
## group_sephealthy.000            -1.01460    0.15575  -6.514
## group_sepno_rejection.000        0.01341    0.10247   0.131
## group_sepacute_rejection.007    -0.10813    0.11873  -0.911
## group_sepno_rejection.007       -0.11783    0.11341  -1.039
## group_sepacute_rejection.090    -0.06124    0.13250  -0.462
## group_sepno_rejection.090        0.15325    0.14585   1.051
## group_sepacute_rejection.365    -0.15761    0.13837  -1.139
## group_sepno_rejection.365        0.08876    0.14741   0.602
## group_sepacute_rejection.biopsy -0.08519    0.13158  -0.647
## 
## Correlation of Fixed Effects:
##              (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## lg10(crtnn)  -0.967                                                
## grp_sph.000  -0.786  0.667                                         
## grp_sp_.000  -0.222  0.008  0.469                                  
## grp_spc_.007 -0.549  0.407  0.607  0.513                           
## grp_spn_.007 -0.608  0.429  0.705  0.688   0.635                   
## grp_spc_.090 -0.678  0.557  0.673  0.463   0.621        0.653      
## grp_spn_.090 -0.691  0.559  0.698  0.537   0.586        0.721      
## grp_spc_.365 -0.670  0.554  0.659  0.444   0.604        0.635      
## grp_spn_.365 -0.702  0.572  0.704  0.531   0.587        0.721      
## grp_spct_r.  -0.484  0.356  0.540  0.463   0.541        0.569      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## lg10(crtnn)                                                     
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.633                                             
## grp_spc_.365  0.649        0.619                                
## grp_spn_.365  0.637        0.715        0.623                   
## grp_spct_r.   0.555        0.522        0.536        0.524

## eGFR, simpler model (biopsy and healthy people not included)
model_data <- data_long %>% 
  filter(!is.na(il1_a_value)) %>% 
  filter(group != 'healthy',
         time_point != 'biopsy')

model_il1a_value_adjEGFR_simple <-  lmer(log10(il1_a_value) ~ log10(egfr) +
                       group*time_point+
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il1a_value_adjEGFR_simple , type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_a_value)
##                       F Df Df.res  Pr(>F)  
## log10(egfr)      0.5022  1 384.49 0.47897  
## group            1.7232  1 124.74 0.19169  
## time_point       2.4330  3 342.57 0.06483 .
## group:time_point 1.5602  3 306.31 0.19906  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR_simple)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_a_value) ~ log10(egfr) + group * time_point + (1 |      id)
##    Data: model_data
## 
## REML criterion at convergence: 692.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5624 -0.2597  0.1086  0.5747  2.0997 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.05608  0.2368  
##  Residual             0.24028  0.4902  
## Number of obs: 424, groups:  id, 138
## 
## Fixed effects:
##                                    Estimate Std. Error t value
## (Intercept)                         0.62772    0.11230   5.590
## log10(egfr)                         0.07336    0.10298   0.712
## groupacute_rejection               -0.01168    0.10338  -0.113
## time_point007                      -0.14357    0.08816  -1.628
## time_point090                       0.11905    0.12944   0.920
## time_point365                       0.05412    0.13118   0.413
## groupacute_rejection:time_point007  0.01947    0.13016   0.150
## groupacute_rejection:time_point090 -0.20129    0.14678  -1.371
## groupacute_rejection:time_point365 -0.23553    0.15074  -1.562
## 
## Correlation of Fixed Effects:
##             (Intr) lg10() grpct_ tm_007 tm_090 tm_365 g_:_00 g_:_09
## log10(egfr)  0.866                                                 
## grpct_rjctn -0.279 -0.008                                          
## time_pnt007 -0.765 -0.585  0.285                                   
## time_pnt090 -0.754 -0.667  0.197  0.614                            
## time_pnt365 -0.761 -0.678  0.194  0.618  0.637                     
## grpct_:_007  0.172 -0.003 -0.654 -0.444 -0.150 -0.148              
## grpct_:_090  0.201  0.053 -0.582 -0.229 -0.524 -0.199  0.461       
## grpct_:_365  0.193  0.049 -0.568 -0.221 -0.193 -0.503  0.449  0.428

There is sign that renal function affect also IL-1b, but it seems to be driven mainly by the difference of healthy subject. Analysis taking into account only planned measurements from patients does not show any indication that renal function affects the IL-1b levels.

4.5.2 IL-1b

Open code
model_data <- data_long %>% 
  filter(!is.na(il1_b_value))

## eGFR
model_il1a_value_adjEGFR <-  lmer(log10(il1_b_value) ~ log2(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il1a_value_adjEGFR, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_b_value)
##                 F Df Df.res  Pr(>F)    
## log2(egfr)  4.287  1 456.46 0.03897 *  
## group_sep  61.156  9 423.86 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_b_value) ~ log2(egfr) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: -505.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.2586 -0.4814 -0.0821  0.2832  7.0701 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.003987 0.06314 
##  Residual             0.015821 0.12578 
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      1.21000    0.03282  36.868
## log2(egfr)                       0.01602    0.00770   2.081
## group_sephealthy.000            -0.65544    0.04072 -16.094
## group_sepno_rejection.000        0.04783    0.02671   1.791
## group_sepacute_rejection.007    -0.07725    0.03074  -2.513
## group_sepno_rejection.007       -0.03984    0.02954  -1.349
## group_sepacute_rejection.090    -0.08160    0.03444  -2.369
## group_sepno_rejection.090       -0.11988    0.03798  -3.157
## group_sepacute_rejection.365    -0.10136    0.03589  -2.825
## group_sepno_rejection.365       -0.13717    0.03834  -3.577
## group_sepacute_rejection.biopsy -0.07194    0.03400  -2.116
## 
## Correlation of Fixed Effects:
##              (Intr) lg2(g) g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log2(egfr)    0.730                                                
## grp_sph.000  -0.865 -0.669                                         
## grp_sp_.000  -0.568  0.007  0.457                                  
## grp_spc_.007 -0.714 -0.420  0.609  0.497                           
## grp_spn_.007 -0.831 -0.427  0.704  0.688   0.632                   
## grp_spc_.090 -0.780 -0.569  0.675  0.444   0.627        0.649      
## grp_spn_.090 -0.812 -0.560  0.700  0.533   0.587        0.725      
## grp_spc_.365 -0.762 -0.563  0.660  0.428   0.610        0.631      
## grp_spn_.365 -0.817 -0.571  0.704  0.528   0.588        0.725      
## grp_spct_r.  -0.634 -0.364  0.541  0.451   0.544        0.565      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log2(egfr)                                                      
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.634                                             
## grp_spc_.365  0.655        0.619                                
## grp_spn_.365  0.637        0.721        0.622                   
## grp_spct_r.   0.559        0.522        0.539        0.523

## creatinine
model_il1a_value_adjCRE <-  lmer(log10(il1_b_value) ~ log10(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il1a_value_adjCRE, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_b_value)
##                         F Df Df.res Pr(>F)    
## log10(creatinine)  2.9475  1 446.93 0.0867 .  
## group_sep         60.3664  9 423.90 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjCRE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_b_value) ~ log10(creatinine) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: -506.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.2579 -0.4712 -0.0827  0.2668  7.0904 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.004023 0.06343 
##  Residual             0.015848 0.12589 
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      1.30700    0.08805  14.843
## log10(creatinine)               -0.05210    0.03020  -1.725
## group_sephealthy.000            -0.64554    0.04068 -15.868
## group_sepno_rejection.000        0.04709    0.02675   1.760
## group_sepacute_rejection.007    -0.07221    0.03066  -2.355
## group_sepno_rejection.007       -0.03548    0.02961  -1.198
## group_sepacute_rejection.090    -0.07411    0.03429  -2.161
## group_sepno_rejection.090       -0.11234    0.03802  -2.955
## group_sepacute_rejection.365    -0.09383    0.03580  -2.621
## group_sepno_rejection.365       -0.12953    0.03843  -3.371
## group_sepacute_rejection.biopsy -0.06735    0.03399  -1.981
## 
## Correlation of Fixed Effects:
##              (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## lg10(crtnn)  -0.967                                                
## grp_sph.000  -0.786  0.667                                         
## grp_sp_.000  -0.222  0.008  0.468                                  
## grp_spc_.007 -0.551  0.412  0.605  0.505                           
## grp_spn_.007 -0.608  0.429  0.705  0.694   0.630                   
## grp_spc_.090 -0.681  0.563  0.672  0.455   0.624        0.648      
## grp_spn_.090 -0.691  0.559  0.699  0.542   0.584        0.726      
## grp_spc_.365 -0.672  0.559  0.658  0.437   0.607        0.631      
## grp_spn_.365 -0.702  0.572  0.704  0.536   0.585        0.726      
## grp_spct_r.  -0.487  0.361  0.539  0.456   0.542        0.564      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## lg10(crtnn)                                                     
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.632                                             
## grp_spc_.365  0.652        0.617                                
## grp_spn_.365  0.635        0.722        0.621                   
## grp_spct_r.   0.557        0.521        0.538        0.522

4.5.3 IL-1a

Open code
model_data <- data_long %>% 
  filter(!is.na(il1_ra_value))

## eGFR
model_il1a_value_adjEGFR <-  lmer(log10(il1_ra_value) ~ log10(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il1a_value_adjEGFR, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_ra_value)
##                  F Df Df.res    Pr(>F)    
## log10(egfr) 0.3599  1 470.56    0.5489    
## group_sep   8.3635  9 421.80 1.647e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_ra_value) ~ log10(egfr) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 328.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.3925 -0.6091 -0.1129  0.4015  3.4088 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.02745  0.1657  
##  Residual             0.08422  0.2902  
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      3.14221    0.07774  40.419
## log10(egfr)                     -0.03650    0.06055  -0.603
## group_sephealthy.000            -0.54051    0.09652  -5.600
## group_sepno_rejection.000        0.11687    0.06332   1.846
## group_sepacute_rejection.007    -0.10337    0.07129  -1.450
## group_sepno_rejection.007       -0.04941    0.07003  -0.706
## group_sepacute_rejection.090    -0.08403    0.08022  -1.047
## group_sepno_rejection.090       -0.03606    0.08969  -0.402
## group_sepacute_rejection.365    -0.25929    0.08351  -3.105
## group_sepno_rejection.365       -0.14709    0.09055  -1.624
## group_sepacute_rejection.biopsy -0.15415    0.07893  -1.953
## 
## Correlation of Fixed Effects:
##              (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log10(egfr)   0.730                                                
## grp_sph.000  -0.864 -0.669                                         
## grp_sp_.000  -0.568  0.007  0.457                                  
## grp_spc_.007 -0.708 -0.429  0.605  0.482                           
## grp_spn_.007 -0.830 -0.427  0.703  0.699   0.622                   
## grp_spc_.090 -0.776 -0.579  0.671  0.429   0.631        0.639      
## grp_spn_.090 -0.814 -0.560  0.700  0.544   0.583        0.735      
## grp_spc_.365 -0.757 -0.571  0.656  0.414   0.614        0.622      
## grp_spn_.365 -0.818 -0.571  0.705  0.539   0.584        0.735      
## grp_spct_r.  -0.630 -0.373  0.538  0.437   0.548        0.556      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log10(egfr)                                                     
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.630                                             
## grp_spc_.365  0.660        0.615                                
## grp_spn_.365  0.634        0.734        0.619                   
## grp_spct_r.   0.563        0.519        0.543        0.520

## creatinine
model_il1a_value_adjCRE <-  lmer(log10(il1_ra_value) ~ log10(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il1a_value_adjCRE, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il1_ra_value)
##                        F Df Df.res    Pr(>F)    
## log10(creatinine) 0.0000  1 462.62    0.9966    
## group_sep         9.0465  9 421.81 1.543e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjCRE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il1_ra_value) ~ log10(creatinine) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 328.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.3829 -0.5912 -0.1104  0.4121  3.4478 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.02792  0.1671  
##  Residual             0.08401  0.2898  
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                      3.1772369  0.2085549  15.235
## log10(creatinine)               -0.0003088  0.0715359  -0.004
## group_sephealthy.000            -0.5796413  0.0964005  -6.013
## group_sepno_rejection.000        0.1171901  0.0633890   1.849
## group_sepacute_rejection.007    -0.1218307  0.0709551  -1.717
## group_sepno_rejection.007       -0.0675076  0.0701644  -0.962
## group_sepacute_rejection.090    -0.1121288  0.0797432  -1.406
## group_sepno_rejection.090       -0.0665761  0.0897262  -0.742
## group_sepacute_rejection.365    -0.2881525  0.0831712  -3.465
## group_sepno_rejection.365       -0.1785399  0.0906953  -1.969
## group_sepacute_rejection.biopsy -0.1720801  0.0787602  -2.185
## 
## Correlation of Fixed Effects:
##              (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## lg10(crtnn)  -0.967                                                
## grp_sph.000  -0.786  0.667                                         
## grp_sp_.000  -0.221  0.008  0.468                                  
## grp_spc_.007 -0.556  0.422  0.601  0.489                           
## grp_spn_.007 -0.608  0.429  0.704  0.706   0.620                   
## grp_spc_.090 -0.686  0.573  0.668  0.439   0.628        0.638      
## grp_spn_.090 -0.692  0.560  0.700  0.554   0.579        0.736      
## grp_spc_.365 -0.676  0.568  0.654  0.422   0.611        0.621      
## grp_spn_.365 -0.703  0.572  0.705  0.548   0.581        0.736      
## grp_spct_r.  -0.491  0.370  0.535  0.442   0.546        0.555      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## lg10(crtnn)                                                     
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.628                                             
## grp_spc_.365  0.657        0.613                                
## grp_spn_.365  0.632        0.735        0.617                   
## grp_spct_r.   0.561        0.517        0.541        0.519

4.5.4 IL-36b

Open code
model_data <- data_long %>% 
  filter(!is.na(il36_b_value))

## eGFR
model_il1a_value_adjEGFR <-  lmer(log10(il36_b_value) ~ log10(egfr) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il1a_value_adjEGFR, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il36_b_value)
##                  F Df Df.res    Pr(>F)    
## log10(egfr) 1.3490  1 483.16     0.246    
## group_sep   8.4876  9 419.12 1.086e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjEGFR)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il36_b_value) ~ log10(egfr) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 568.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3808 -0.3381  0.0098  0.4113  3.1354 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.05718  0.2391  
##  Residual             0.13051  0.3613  
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      0.79421    0.10017   7.929
## log10(egfr)                      0.09076    0.07778   1.167
## group_sephealthy.000            -0.61739    0.12452  -4.958
## group_sepno_rejection.000        0.05236    0.08192   0.639
## group_sepacute_rejection.007    -0.21499    0.08931  -2.407
## group_sepno_rejection.007       -0.09652    0.09049  -1.067
## group_sepacute_rejection.090    -0.43762    0.10105  -4.331
## group_sepno_rejection.090       -0.21276    0.11517  -1.847
## group_sepacute_rejection.365    -0.60102    0.10508  -5.720
## group_sepno_rejection.365       -0.26771    0.11628  -2.302
## group_sepacute_rejection.biopsy -0.55169    0.09899  -5.573
## 
## Correlation of Fixed Effects:
##              (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## log10(egfr)   0.728                                                
## grp_sph.000  -0.863 -0.666                                         
## grp_sp_.000  -0.570  0.006  0.458                                  
## grp_spc_.007 -0.700 -0.439  0.599  0.462                           
## grp_spn_.007 -0.829 -0.425  0.701  0.715   0.608                   
## grp_spc_.090 -0.768 -0.591  0.665  0.409   0.637        0.625      
## grp_spn_.090 -0.815 -0.559  0.701  0.561   0.576        0.748      
## grp_spc_.365 -0.750 -0.581  0.650  0.395   0.620        0.608      
## grp_spn_.365 -0.820 -0.570  0.705  0.555   0.578        0.748      
## grp_spct_r.  -0.624 -0.384  0.532  0.418   0.552        0.544      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## log10(egfr)                                                     
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.624                                             
## grp_spc_.365  0.666        0.609                                
## grp_spn_.365  0.628        0.750        0.613                   
## grp_spct_r.   0.568        0.514        0.548        0.515

## creatinine
model_il1a_value_adjCRE <-  lmer(log10(il36_b_value) ~ log10(creatinine) +
                       group_sep +
                       (1|id), data = model_data, REML = TRUE)

Anova(model_il1a_value_adjCRE, type=2, test.statistic='F')
## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: log10(il36_b_value)
##                        F Df Df.res    Pr(>F)    
## log10(creatinine) 1.4245  1 476.83    0.2333    
## group_sep         8.5998  9 419.39 7.347e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_il1a_value_adjCRE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log10(il36_b_value) ~ log10(creatinine) + group_sep + (1 | id)
##    Data: model_data
## 
## REML criterion at convergence: 567.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3756 -0.3412  0.0117  0.4202  3.1371 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 0.05725  0.2393  
##  Residual             0.13044  0.3612  
## Number of obs: 500, groups:  id, 186
## 
## Fixed effects:
##                                 Estimate Std. Error t value
## (Intercept)                      1.01992    0.26817   3.803
## log10(creatinine)               -0.11026    0.09196  -1.199
## group_sephealthy.000            -0.61969    0.12431  -4.985
## group_sepno_rejection.000        0.05102    0.08192   0.623
## group_sepacute_rejection.007    -0.21531    0.08895  -2.421
## group_sepno_rejection.007       -0.09802    0.09058  -1.082
## group_sepacute_rejection.090    -0.43842    0.10049  -4.363
## group_sepno_rejection.090       -0.21481    0.11516  -1.865
## group_sepacute_rejection.365    -0.60224    0.10471  -5.752
## group_sepno_rejection.365       -0.27011    0.11640  -2.320
## group_sepacute_rejection.biopsy -0.55244    0.09882  -5.590
## 
## Correlation of Fixed Effects:
##              (Intr) lg10() g_.000 g__.000 grp_spc_.007 grp_spn_.007
## lg10(crtnn)  -0.967                                                
## grp_sph.000  -0.784  0.665                                         
## grp_sp_.000  -0.222  0.008  0.468                                  
## grp_spc_.007 -0.561  0.432  0.595  0.470                           
## grp_spn_.007 -0.607  0.427  0.702  0.720   0.606                   
## grp_spc_.090 -0.692  0.585  0.663  0.420   0.634        0.625      
## grp_spn_.090 -0.693  0.559  0.701  0.568   0.573        0.749      
## grp_spc_.365 -0.681  0.578  0.648  0.404   0.616        0.608      
## grp_spn_.365 -0.703  0.571  0.706  0.562   0.575        0.749      
## grp_spct_r.  -0.497  0.381  0.531  0.424   0.550        0.543      
##              grp_spc_.090 grp_spn_.090 grp_spc_.365 grp_spn_.365
## lg10(crtnn)                                                     
## grp_sph.000                                                     
## grp_sp_.000                                                     
## grp_spc_.007                                                    
## grp_spn_.007                                                    
## grp_spc_.090                                                    
## grp_spn_.090  0.622                                             
## grp_spc_.365  0.664        0.607                                
## grp_spn_.365  0.627        0.750        0.612                   
## grp_spct_r.   0.566        0.512        0.546        0.514

5 IL18 forms plot

At first, we need to transform IL-18 levels to molar concentration (pM/L)

Open code
data_long <- data_long %>% 
  mutate( il18bp_pM = (((il18_bp_value/1e12)*1000) / 19000)*1e12,
         il18tot_pM = (((il18_value/1e12)*1000) / 17200)*1e12,
        il18free_pM = (((il18_free/1e12)*1000) / 17200) *1e12)

Now we can construct plots showing the relationship between IL-18BP levels and levels of IL-18 (free and total). We will show it with focus on healthy individuals, and for 1st measurement of the group without rejection where the IL-18BP and IL-18 were most variable, with examples of minimal as well as vy high bounding

Open code
## prepare data
il18d <- data_long %>% 
  filter(!is.na(il18_value),
         !is.na(il18_bp_value),
         !is.na(il18_free)) %>% 
  mutate(group = factor(group, levels = 
                          c("no_rejection", 
                            "acute_rejection",
                            "healthy")),
         time_point = factor(time_point)) %>% 
  mutate(time_point = if_else(group == 'healthy', 'healthy', time_point)) %>% 
  mutate(group_sep = interaction(group, time_point))

## Plot for healthy
p18_healthy <- il18d %>% 
  filter(group_sep == 'healthy.healthy') %>% 
  ggplot(aes(x = il18bp_pM)) +
  geom_point(aes(y = il18free_pM, color = "IL-18 Free")) +  
  geom_point(aes(y = il18tot_pM, color = "IL-18 Total"), pch=1) +  
  geom_line(data = . %>% 
              select(il18bp_pM, il18free_pM,  il18tot_pM) %>% 
              pivot_longer(cols = -il18bp_pM, 
                           names_to = "variable", values_to = "value"),
            aes(x = il18bp_pM, y = value, group = il18bp_pM), col = 'grey20') +
  scale_color_manual(values = c("IL-18 Free" = "blue", "IL-18 Total" = "red")) +
  labs(x = "IL-18BP (pM)", y = "IL-18 (pM)", color = "") +
  coord_cartesian(ylim = c(0, NA)) +
  ggtitle("Samples from healthy donors")
p18_healthy

Open code
if(file.exists('gitignore/figures/p18_healthy.pdf') == FALSE){
  ggsave("gitignore/figures/p18_healthy.pdf", 
         plot = p18_healthy, width = 8, height = 5.5, units = "in")
}

## plot for D0, group without rejection
p18_no_reject_0 <- il18d %>% 
  filter(group_sep == 'no_rejection.000') %>% 
  ggplot(aes(x = il18bp_pM)) +
  geom_point(aes(y = il18free_pM, color = "IL-18 Free")) +  
  geom_point(aes(y = il18tot_pM, color = "IL-18 Total"), pch=1) +  
  geom_line(data = . %>% 
              select(il18bp_pM, il18free_pM,  il18tot_pM) %>% 
              pivot_longer(cols = -il18bp_pM, 
                           names_to = "variable", values_to = "value"),
            aes(x = il18bp_pM, y = value, group = il18bp_pM), col = 'grey20') +
  scale_color_manual(values = c("IL-18 Free" = "blue", "IL-18 Total" = "red")) +
  labs(x = "IL-18BP (pM)", y = "IL-18 (pM)", color = "") +
  coord_cartesian(ylim = c(0, NA)) +
  ggtitle("No rejection group, the 1st measurement")
p18_no_reject_0

Open code
if(file.exists('gitignore/figures/p18_no_reject_0.pdf') == FALSE){
  ggsave("gitignore/figures/p18_no_reject_0.pdf", 
         plot = p18_no_reject_0, width = 8, height = 5.5, units = "in")
}

## plot for D0, group with rejection
p18_acute_reject_0 <- il18d %>% 
  filter(group_sep == 'acute_rejection.000') %>% 
  ggplot(aes(x = il18bp_pM)) +
  geom_point(aes(y = il18free_pM, color = "IL-18 Free")) +  
  geom_point(aes(y = il18tot_pM, color = "IL-18 Total"), pch=1) +  
  geom_line(data = . %>% 
              select(il18bp_pM, il18free_pM,  il18tot_pM) %>% 
              pivot_longer(cols = -il18bp_pM, 
                           names_to = "variable", values_to = "value"),
            aes(x = il18bp_pM, y = value, group = il18bp_pM), col = 'grey20') +
  scale_color_manual(values = c("IL-18 Free" = "blue", "IL-18 Total" = "red")) +
  labs(x = "IL-18BP (pM)", y = "IL-18 (pM)", color = "") +
  coord_cartesian(ylim = c(0, NA)) +
  ggtitle("Acute rejection group, the 1st measurement")
p18_acute_reject_0

Open code
if(file.exists('gitignore/figures/p18_acute_reject_0.pdf') == FALSE){
  ggsave("gitignore/figures/p18_acute_reject_0.pdf", 
         plot = p18_acute_reject_0, width = 8, height = 5.5, units = "in")
}

6 Reproducibility

Open code
sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=cs_CZ.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=cs_CZ.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=cs_CZ.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=cs_CZ.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Prague
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lqmm_1.5.8       effects_4.2-2    rstudioapi_0.13  quantreg_5.88   
##  [5] SparseM_1.81     coxed_0.3.3      survival_3.5-7   rms_6.7-1       
##  [9] Hmisc_5.1-1      bayesplot_1.8.1  ggdist_3.3.0     kableExtra_1.3.4
## [13] lubridate_1.8.0  corrplot_0.92    arm_1.13-1       lme4_1.1-34     
## [17] MASS_7.3-60      janitor_2.2.0    RJDBC_0.2-10     rJava_1.0-6     
## [21] DBI_1.1.2        projpred_2.6.0   brms_2.16.3      Rcpp_1.0.11     
## [25] glmnet_4.1-8     Matrix_1.6-3     boot_1.3-28      cowplot_1.1.1   
## [29] pROC_1.18.0      mgcv_1.9-1       nlme_3.1-163     openxlsx_4.2.5  
## [33] flextable_0.9.3  sjPlot_2.8.15    car_3.1-2        carData_3.0-5   
## [37] gtsummary_1.7.2  emmeans_1.7.2    ggpubr_0.4.0     forcats_1.0.0   
## [41] stringr_1.5.0    dplyr_1.1.3      purrr_1.0.2      readr_2.1.2     
## [45] tidyr_1.2.0      tibble_3.2.1     ggplot2_3.4.3    tidyverse_1.3.1 
## 
## loaded via a namespace (and not attached):
##   [1] fs_1.6.3                matrixStats_1.0.0       httr_1.4.2             
##   [4] webshot_0.5.5           insight_0.19.5          SparseGrid_0.8.2       
##   [7] tools_4.3.2             backports_1.4.1         sjlabelled_1.2.0       
##  [10] utf8_1.2.3              R6_2.5.1                DT_0.17                
##  [13] withr_2.5.0             Brobdingnag_1.2-7       prettyunits_1.1.1      
##  [16] gridExtra_2.3           cli_3.6.1               textshaping_0.3.6      
##  [19] performance_0.10.5      gt_0.9.0                shinyjs_2.1.0          
##  [22] officer_0.6.2           sandwich_3.0-1          labeling_0.4.2         
##  [25] sass_0.4.7              mvtnorm_1.1-3           polspline_1.1.24       
##  [28] ggridges_0.5.3          askpass_1.1             commonmark_1.9.0       
##  [31] systemfonts_1.0.4       StanHeaders_2.21.0-7    foreign_0.8-86         
##  [34] gfonts_0.2.0            svglite_2.1.0           readxl_1.3.1           
##  [37] httpcode_0.3.0          generics_0.1.2          shape_1.4.6            
##  [40] gtools_3.9.2            crosstalk_1.2.0         distributional_0.3.2   
##  [43] zip_2.2.0               inline_0.3.19           loo_2.4.1              
##  [46] fansi_1.0.4             abind_1.4-5             lifecycle_1.0.3        
##  [49] multcomp_1.4-18         yaml_2.3.5              snakecase_0.11.1       
##  [52] grid_4.3.2              promises_1.2.0.1        crayon_1.5.2           
##  [55] miniUI_0.1.1.1          lattice_0.22-5          haven_2.4.3            
##  [58] pillar_1.9.0            knitr_1.37              estimability_1.3       
##  [61] shinystan_2.5.0         codetools_0.2-19        glue_1.6.2             
##  [64] fontLiberation_0.1.0    data.table_1.14.8       broom.helpers_1.14.0   
##  [67] vctrs_0.6.3             cellranger_1.1.0        gtable_0.3.4           
##  [70] assertthat_0.2.1        xfun_0.40               mime_0.12              
##  [73] survey_4.2-1            coda_0.19-4             rsconnect_0.8.25       
##  [76] iterators_1.0.14        shinythemes_1.2.0       ellipsis_0.3.2         
##  [79] TH.data_1.1-0           pbkrtest_0.5.2          xts_0.12.1             
##  [82] fontquiver_0.2.1        threejs_0.3.3           rstan_2.21.3           
##  [85] tensorA_0.36.2          rpart_4.1.23            colorspace_2.1-0       
##  [88] nnet_7.3-19             tidyselect_1.2.0        processx_3.8.2         
##  [91] compiler_4.3.2          curl_4.3.2              rvest_1.0.2            
##  [94] htmlTable_2.4.0         xml2_1.3.3              fontBitstreamVera_0.1.1
##  [97] bayestestR_0.13.1       colourpicker_1.1.1      posterior_1.2.0        
## [100] checkmate_2.0.0         scales_1.2.1            dygraphs_1.1.1.6       
## [103] callr_3.7.3             digest_0.6.29           minqa_1.2.4            
## [106] rmarkdown_2.25          htmltools_0.5.6         pkgconfig_2.0.3        
## [109] base64enc_0.1-3         highr_0.9               dbplyr_2.1.1           
## [112] fastmap_1.1.0           rlang_1.1.1             htmlwidgets_1.6.2      
## [115] shiny_1.5.0             farver_2.1.1            zoo_1.8-9              
## [118] jsonlite_1.7.3          magrittr_2.0.3          Formula_1.2-4          
## [121] munsell_0.5.0           gdtools_0.3.3           stringi_1.7.12         
## [124] plyr_1.8.8              pkgbuild_1.3.1          parallel_4.3.2         
## [127] sjmisc_2.8.9            ggeffects_1.3.1         splines_4.3.2          
## [130] hms_1.1.1               sjstats_0.18.2          ps_1.6.0               
## [133] igraph_1.5.1            uuid_1.0-3              ggsignif_0.6.3         
## [136] markdown_1.8            reshape2_1.4.4          stats4_4.3.2           
## [139] rstantools_2.1.1        crul_1.4.0              reprex_2.0.1           
## [142] evaluate_0.15           mitools_2.4             RcppParallel_5.1.7     
## [145] modelr_0.1.8            nloptr_2.0.0            tzdb_0.4.0             
## [148] foreach_1.5.2           httpuv_1.6.5            MatrixModels_0.5-0     
## [151] openssl_1.4.6           broom_1.0.5             xtable_1.8-4           
## [154] rstatix_0.7.0           later_1.3.0             viridisLite_0.4.2      
## [157] ragg_1.2.1              cluster_2.1.6           gamm4_0.2-6            
## [160] bridgesampling_1.1-2